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1. Introduction. Through a consideration of the fundamental aspects of the uni- 
verse Newton was led to the invention of fluxions. The physical idea of “rate of 
change,” the geometric idea of “slope of a curve” together with his newly invented 
mathematics proved to be very powerful in describing and predicting a wide range 
of phenomena of nature. Since Newton's day, an enormous number of physical phe- 
nomena have been described in terms of a few “laws of nature.” Very often these 
laws make use of the calculus, especially when applied to a specific problem. Thus 
large sections of the phenomena of the physical universe are described by the solutions 
of differential equations for the appropriate boundary conditions. 

The engineer in his attempt to make nature work his way is continually presented 
with problems to be solved. These problems, even when very technical in nature, 
usually contain an element not present in problems considered by mathematicians or 
physicists. A “solution” of an engineer's problem is often a numerical answer or a 
graph obtained in a specified time. A poor answer which meets the deadline date is 
far superior to a precise answer a week later. Thus the engineer, or any applied sci- 
entist, should, when choosing the method of attack on a problem, keep in mind the 
time when the answer is due. 

The following methods of solution are available: 

1) the answer can be guessed; 

2) some experiments can be run; 

3) the result can be computed from the basic laws of nature, by use of whatever 
mathematical methods are needed. 

Obviously the first method has one certain and everlasting superiority over the 
second and third; it is quick. It can meet any deadline set at a future time. It, of 
course, has one big disadvantage. The solution is always of doubtful quantitative 
value even though it is of immense qualitative value. In fact, the guessing process 
should always be used as a guide to correct results by methods 2 or 3. 

If time permits, the answer may be sought by experiment or computation, or both. 
At the present time computational methods, when applied to real physical systems 
to obtain solutions of sufficient accuracy, are often so cumbersome that the vast 
majority of engineering problems are solved primarily by experiment. Usually com- 
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putations are confined to one dimensional approximations to the real problem, or to 
certain simple two or three dimensional approximations. 

In the field of partial differential equations, which includes the key to a vast 
number of practical problems, the present day need is for methods of solution speedier 
and more general than the majority of analytical methods thus far discovered. Nu- 
merical methods of solution aim at securing an approximate quantitative answer to a 
given problem as directly as possible from the statement of the problem, so as to 
reduce the time required for solution. An advantage of this more direct approach is the 
possibility of injecting into the solution any facts known (or supposed known) by the 
computer. Thus, physical “facts” of any nature can be made use of, and if some of the 
“facts” are wrong the solution will indicate this. Obviously, this has an enormous 
advantage over analytical methods where generally little use can be made of detailed 
physical observations, except as a check on the final result. 

2. An elementary problem. By way of illustrating the methods involved, let us 
consider the transient flow of heat in a two-dimensional homogeneous, isotropic solid 
for which the physical law of the conservation of energy yields the familiar differential 
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equation 


where J denotes the temperature, x and y space coordinates, ¢ the time, a the ther- 
mal diffusivity (a property of the material, assumed constant). The numerical meth- 
ods considered here are based upon a finite difference approximation to the differential 
equation. The simplest of these may be obtained as follows. 
By definition 
dT (x, ¥, t) T(x, vy, t + 6t) — T(x, y, d 
———— = lim - - (2) 
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and similar expressions hold for the derivatives with respect to x and y. As an ap- 
proximation, the limit operation may be omitted. If this is done the following expres- 


sions are obtained 
0°T T(x — bx, y, t) — 2T(x, y, t) + T(x + dz, y, 2) T3; —27,+ 7; (3) 
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where the last expressions on the right follow the notation of Fig. 1. Substitution of 
these approximations into Eq. (1) gives, after some rearrangement, 








j abt aét 2abt 2aét 
T(x, y, t + 64) = —(T3 + T:) + —(T2 + T,) +(1 — — \re (5) 
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which relates the temperature at a given point 0 at time ‘+ 6¢ to the temperatures 
which existed in the neighborhood of 0 at time ¢. Since the only restriction on 5t, 6x, dy 
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is that they be small enough to render the finite difference approximation sufficiently 
accurate, we set 

4a65t = bx? = by? = 6, (6) 
whence 


T(x, y,t + 6) = 2(T1 + T2 + T3 + Ty). (7) 


This same equation can be derived directly from the physical problem by making 
physical assumptions only. Such an approach is very valuable, since many more or 
less vague physical assumptions are always made before any engineering problem can 
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be set up mathematically. Let us assume that the material of the solid body is divided 
up into three sets of overlapping squares (or rectangles if desired). One set represents 
all heat conduction of the body in the x-direction, as indicated by the dotted square 
between points 0 and 1 of Fig. 1. The secord set of squares represents the conduction 
in the y-direction and finally the third set, surrounding each point, represents all the 
material as far as heat capacity is concerned. The thermal energy Qj~» which is con- 
ducted in unit time to point 0 along the rod 1—0 is obtained from Fourier’s heat 
conduction equation 


A dT (8) 
i dx 
in the form 
Qi-o = kb(T1 — To), (9) 


where b is the thickness of the two dimensional body considered. The energy arriving 
at the point 0 from all the surrounding points is thus 


Qo = kb(T1 + T2 + Ts + Ts — 47). (10) 


This heat will result in an increase of temperature of the material associated with 


point 0, whence 
, 2 (%, y, t + dt) —_ To 
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Qo = cpbé (11) 
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We equate (10) and (11) and rearrange terms. Thus, 





cpd” 


kit 4 kit 
T(x, y, t + 5t) ote ree T3 + ri) +(1- ; )t (12) 
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which reduces to (7), since the thermal diffusivity a =k/cp, and 6t, 6 are related by (6). 

The method of use of Eq. (7) for the solution of a transient heat flow problem 
expressed by the differential Eq. (1) is direct. The space domain is divided into a 
square net of points. (The exact relationship between the required net spacing and the 
desired accuracy of solution has not yet been studied, to the author’s knowledge.) 
The initial values (at ¢=0) of temperature are attached to each point and the values 
at successive times are computed by the averaging process. Figs. 2 illustrates such a 
solution. 
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Fic. 2a. (See legend below Fig. 2c.) 


By the nature of the process it is clear that the shape of the domain and the 
boundary conditions (generally some relation between the boundary value and the 
normal derivative) cause no special difficulty, such as occurs in the analytical ap- 
proach, since they can be treated numerically as required for the points nearest (or on) 
the boundary. Boundary conditions and net spacing should be chosen of comparable 
accuracy as judged physically (in the absence of rigorous methods). Experience indi- 
cates that it is seldom worth the trouble to derive special boundary formulas, since 
generally linear extrapolation or a simple plot on graph paper will give the same ac- 
curacy much more speedily. 

Let us return to Eq. (1) and the attendant physical problem. It is known that if 
the thermal conditions at the physical boundaries are held fixed (in time), the in- 
terior temperature will (after a sufficient period) reach a steady value (to any given 
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degree of approximation). In the present problem, this means that with a sufficiently 
large number of applications of Eq. (7) the solution of Laplace’s equation will be 
obtained for the given domain and boundary conditions. The process proposed above, 
the one closely followed by nature, is numerically very slow unless an elaborate com- 
puting machine is used [1 ].* 

Several schemes have been proposed for speeding up this process. By starting 
with Laplace’s equation and making the finite difference approximations, by starting 
physically with Eq. (10) and Qo=0, or by observing that for steady conditions 
T(x, y, t+6t) = To in Eq. (7), one obtains 


To = 3(71 + T2 + Ts + Ty). (13) 


When Eg. (13) is written for each of m points of a square net covering a given 
domain, m linear equations result. Several iteration processes have been devised for 
solving such a system [2, 5, 6], and the convergence of some of these has been 
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discussed [3, 4, 7, 8]. All of these methods propose computations on the values of T 
by a specifically stated iteration process which can be shown rigorously to converge 
and can be performed by a completely automatic computing machine. In no case 
is it possible to add physical information after the values guessed initially are at- 
tached to each net point, without upsetting the scheme of solution. 

A new scheme for solving Eq. (13) for the m points of a two-dimensional domain 
is given by the relaxation method [9]. This method is so superior to others in point 
of the time required to reach a solution of given accuracy that it will be discussed in 


* The numbers in square brackets refer to the bibliography at the end of the paper. 
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detail. Again, the physical problem will be considered as one of heat conduction, 
although any other phenomena leading to Laplace’s equation would serve as well. 
R. V. Southwell [10], who developed the relaxation method from a consideration of 
problems of statics, generally speaks in terms of a tension net as an approximation 
to a soap film or membrane. 
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Fics. 2a-c. Transient heating of furnace wall (Fig. 5). At time 0 wall at uniform temperature 
T=100° F has inner surface temperature raised to 500° F. Wall thermal diffusivity a=.01 ft?/hr. 
6: =.208 ft. Therefore 5;=1.083 hrs. The 6 temperatures shown at each point are at time intervals of 
1.083 hrs. Each temperature is the average of the temperatures at the 4 surrounding points at the pre- 


ceding time. Heat loss is }_Q’ along inner surface. 
(From The numerical solution of heat conduction problems. By H. W. Emmons, Trans. A.S.M.E., 


65, 607-615 (1943)). 


Instead of focusing attention on the values of JT and the averaging process of 
Eq. (13), let us return to Eq. (10). To solve a problem the domain is drawn and the 
net points chosen. Values of T are then attached (by guess or any information avail- 
able from experiment, special solutions, prior work, field mapping, etc.) to each point. 
From these values the residuals 


Q = Qo/kb = T3 + T2 + T3 + Ts — 470 (14) 


are computed and recorded. The Q, thus computed, can be thought of as interior heat 
sinks which must be removed. Now, instead of setting up a specific iteration process, 
we merely observe that if the temperature at one point (0) is altered, all others re- 
maining fixed, the residuals will change according to the pattern in Fig. 3, the “relaxa- 
tion pattern.” Each change of T, at any point, effects a redistribution of the residuals, 
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Q, among the net points, and such changes of T are desired which will move all the 
sinks to the boundary. 
These “operating instructions” may appear vague. Indeed they 
are vague. Their vagueness is the source of their great power be- oy 
cause the computer may without any effort alter the procedure to 4] 
attain more rapid approach to the final answer (of no residuals). as co 
There is only one way to appreciate fully the meaning of these re- ea 
marks and that is to do a problem. 
j é ‘ Fie. 3. 
Let us consider the extremely simple problem of Fig. 4, where 
the solution of Laplace’s equation is desired with zero boundary values. For illustra- 
tions of the various methods, only three interior points are used and the ridiculous 
tfial solution indicated in Fig. 4 is assumed. In this simple problem all the methods are 
equally easy. The real time-saving advantages 
Sante © 0 © of the relaxation process only appear with more 
net points. In Table IA the transient solution 
[differential equation (1)] is carried out by use 
of the difference equation (7). In Table IB 
Liebman’s method [2] is used, Eq. (13). In both 
cases, each value after the first at any point is 
obtained by adding four numbers and dividing 
by 4. This process cannot (for say 4 digit num- 
bers) be carried out mentally. Fifteen and ten 
changes respectively were needed to make the 
error less than unity. 
Tables IC, D, E show applications of the 
relaxation method and these will be followed in 
Fic. 4. detail to illustrate various “tricks” which serve 
to speed up the elimination of residuals. In 
Tables IC, D, E the initial value of T and the subsequent corrections are shown at 
the right of each field, the values of the residuals (heat sinks) are shown to the left. 
The largest heat sink occurs in the vicinity of the greatest deviation of the assumed 
values from the correct solution, so changes are first made at this point. For purposes 
of illustration, each time any change is made in Tables IC, D, E a space is left at all 
unaffected points, so that the work can be followed. (Generally this is not done.) 
Let us consider Table IC. To eliminate exactly a residual of —420 at point 3 
would require a change in 7; of —105. The making of this change is equivalent ex- 
actly to the averaging process. But why do we bother with three digit numbers? 
Since our residuals are very large, let us make simple large changes of about the right 
size. Therefore, let us change 7; by —100. By the relaxation pattern of Fig. 3 the 
residuals become Q; = —420+4(100) = —20, Q.= —80—100 = —180. Now Q: is larg- 
est. Accordingly, we change 72; a change of —50 was chosen and the third residuals 
at each point computed. We notice that the residual at the point 3, much improved 
at the first step, has been spoiled again by the following change at point 2. This al- 
ways happens when a point is surrounded by other points with residuals of the same 
sign. We would have done better to overshoot zero at point 3 on the first change. On 
the next change, at point 1, we overshoot the zero and make the residual positive. The 
residual at point 2 changes from +20 to —20 and the previous change, might well 
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have been larger. By not overshooting enough or by overshooting too much we do 

no harm except that some time is lost. A little practice improves one’s guessing. 
In Table ID the block relaxation 
TABLE I process is used, i.e., a block of points 
A) Solution by Eq. (17); 15 changes. is changed simultaneously. This is de- 
B) Solution by Eq. (13) (Liebman’s method); 10 sirable since all the residuals are nega- 
changes. ; tive. So large a change is made that their 

C) Illustrates overshooting; 5 changes. (f =overshot, ‘ ; : 
a average resultant residual is about zero 
=should have overshot). : : : > 

(or overshot if desired). A change of 50 


D) Illustrates block relaxation; 6 changes made, 
first 3 at once. was chosen. If one computes by the 





























E) Illustrates use of prior knowledge; 3 changes. relaxation pattern [Eq. (14)] directly 
—— ee oe gain by block relaxatiorf. 
Point | i 2 3 | Instead let us consider the physical ar- 
| 0 ~ | 490 a rangement. If the temperatures at two 
15 40 15 adjacent points are both changed by the 
A | 10 7.5 | 10 same amount there will be no change 
1.4 5 1.4 of heat flow along the connecting rod. 
1.2 7 oe Thus the relaxation pattern —1, 1 as 
: - * ____**_| given by Eq. (9) is used for each rod 
40 60 120 independently, and at the points of 
15 33.8 8.4 Table ID, the residuals change by Q;= 
B | 8.4 4.2 1.0 —100+3(50) =50, Q2 = —80+2(50) =20, 
| 1.0 S| ‘1 |) Q3= —420+3(50) = —270. From this 
_| i" a point the solution is continued as before. 
|— 100 40'—- 80 60\|—420 120 Of course we know the solution in 
; — —180 — 20 —100 the present case is zero everywhere. This 
C |-150 20 —50 |— 70 knowledge is used in Table IE immedi- 
10 —404/— 20 , ot ately and the solution obtained in three 
‘* = 7 ~10 7 = steps. Naturally one would not have 
e __| started this problem with so poor an 
1-100 40/— 80 60|—420 120 assumed set of initial values. Also the 
50 —5S0 20-50 |—270 — S0| prior knowledge is never as extensive as 
- 4 = > | =a in the present trivial problem, but ob- 
40 — 10 —10 0 ‘ ‘ : 
0 +10 0 va | servation of trends in the solution often 
gives a clue to the way in which values 
i-100 40/— 80  60/-—420 120] should be changed. 
, — 200 60 —120 It is not possible to judge the speed 
| 160 ao —0 0 | and ease of the relaxation method com- 
0 —40 0 _— . 
Wize _——_ se | pared to the averaging methods by 


comparison of the various procedures 
illustrated in Table I. The relaxation method is far superior when a large number of 
points are used because of its flexibility (possibility of overshooting), its use of simple 
large numbers when the residuals are large and simple small numbers when they have 
become small, and the fact that the most difficult operation is to multiply a simple 
number by 4 and add the result to or subtract it from another, all of which can easily 


be done mentally. 
Figs. 5 illustrate the solution of a somewhat more difficult problem. It should be 
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noted that the solution is started with a few points only. Then, for greater accuracy 
more points are added. These need not be added everywhere but only where greater 
accuracy is desired, or where some variable changes rapidly. In the illustration the 
points were doubled in number by adding one point at the center of each four points. 
The resultant net is square but diagonal to the first. Since Laplace’s equation is in- 
variant on rotation of the axes, the same average formula applies to values on the new 
net. Guessed values for starting the finer net are obtained by averaging the surround- 
ing four values. If a fine net is used locally, as might have been done near the reentrant 
corner in the illustration, it can always be connected with the coarser net through use 
of the diagonal formula where the two nets meet. 
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Fic. 5a. Furnace, showing section under consideration. 


The results of computation of the furnace wall conduction problem are sum- 
marized in Table II. We note particularly the time required for a solution and the 


TABLE II 


Numerical Solution of Heat Conduction 
Problem of Fig. 5 











Number of Calculated Deviation from Hours 
points thermal experimental required for 
used resistance solution calculation 
para 0735 a | 
by arithmetic mean area — | 15.3% -0S 
2b | 
.0806 
12 —— 5 .75 
kb % 
.0825 
19 = 2.8% | 1.75 
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Fic. 5b. Start solution with linear temperature distribution, relax Q’ to get final temperature. Heat 
transferred is Q=kb[232 +208+(4.5/5)202 |] =622 kb. Thermal resistance R=AT/8Q =400/(8 X622kb) 


= .0806/kb. Superscripts indicate step of calculation. °=original values, !:*¢te- =successive steps. 
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Fic. 5c. Start with solution of coarse net Fig. 3a. Relax Q’ to get final temperature. Heat transferred 
Q=kb[195/2+2111+2102+(4/5)102]=606 kb. Thermal Resistance R=AT/8Q =400/(8 X606kb) 
= .0825/kb. 

(From The numerical solution of heat conduction problems. By H. W. Emmons, Trans. A.S.M.E., 
65, 607-615 (1943)). 
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corresponding error. This same problem took 11 hours when solved by the averaging 


process. 

Another problem, still more difficult, is illustrated in Fig. 6. This is the problem 
of the water tube of a boiler, and the boundary conditions are varied. The water in- 
side the tube is assumed to maintain a constant wall temperature T=0. The lower 
half of the tube is embedded in insulating material, so the normal gradient of tempera- 
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Fic. 6. Heat transfer by radiation through furnace water wall tube. 
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ture is there assumed to be zero. On the upper half of the tube heat is transferred by 
radiation in such a way that the heat input is constant per unit projected area. The 
normal gradient of temperature is thus proportional to sin 8. This problem was easiest 
to solve by a transformation from the x, y-plane to the z, 6-plane, where 6=tan-'y/x, 
z=log r=} log (x?+y?). The net points used are indicated in the x, y-plane by the in- 
tersection of the radial lines and circles. The numerical solution checked “exactly” 
with the analytical solution, as observed by superposing the graphed results. From 
an engineering point of view this solution is “exact,” since it deviates much less from 
the analytical solution than the uncertainty of the boundary conditions. 
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So far we have dealt at length with the relation between the well known averaging 
process and the relaxation process. Indeed, it may appear that we have wasted time 
with trivialities of arithmetic but the author’s conversations with many have indi- 
cated that it is just these minor details in point of view that make the relaxation 
process about five times as rapid as the iteration processes. To appreciate fully the 
power of the flexibility of the relaxation method, one must take pencil and paper and 
carry out the numerical process in all its uninteresting details. In fact, for the comp- 
puter (as opposed to those who think only about the logic behind the computation 
methods) the relaxation method has a spirit lacking entirely from the iteration proc- 
esses. The former challenges one’s intellect at each step to make the best possible 
guess, while the latter reduces one to the status of an automatic computing machine 
(without the advantage of no computational errors). It should not be inferred that the 
relaxation process requires high intellectual powers. If changes are chosen in a specifi- 
able way it reduces exactly to the iteration process. The computer can then vary from 
this completely specified process by whatever amount fits his own skill. 

3. Other types of equations solved by the relaxation method. After the essential 
idea of the relaxation method is grasped, other problems may be solved by rather ob- 
vious steps. The differential equation is converted into a difference equation, some 
quantity (to be zero for the solution) is chosen as the residual, and the relaxation pat- 
tern is set up. Changes which make the residuals smaller are then made. We notice 
that no question of convergence can occur in the general relaxation process, since no 
specific instructions are given. If, after some steps, the residuals get worse, the intelli- 
gent computer goes back and makes changes in the opposite direction. These remarks 
oversimplify the problem somewhat because of two facts; first, the computer may 
become confused as to whether or not the residuals are really better, and secondly 
there is always a question of whether or not a solution with zero residuals exists 
(see [8]). 

The following is a partial list of types of problems solved, with some brief details 
of their solution. 


A. Poisson's equation, 
Ag = Giz t Gy = f(%, y). (15) 
(Subscripts denote partial differentiation.) In finite difference form, for a square 
net of points the approximating difference equation is (see Fig. 1 for notation) 


01 + oo + 3 + ¢4 — 460 — I . # y)bx? a Q, (16) 


where the residual Q is to be zero at each net point. Since f(x, y) has a known numeri- 
cal value at each net point, it enters into the value of the residuals at the start of the 
computation, and in no way affects the relaxation pattern used for Laplace’s equation, 
Fig. 3. Thus the solution of Poisson’s equation is precisely as easy as that of Laplace’s. 
B. Biharmonic equation, 
O*w Ow O*w 
AAw = + 2——_- + — = 0. (17) 


2 : 
Ox' 0°xd*y = ay" 





The finite difference equation for a square net of points is most simply derived by 


converting the two A forms separately. Thus 


Aw; + Awe + Aw; + Aw, — 4Auo = 0 (18) 
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from which, by expansion of the remaining A’s, one obtains the relaxation pattern 


shown on Fig. 7, where the solution of a problem of the deflection of a plate with 
clamped deflected edges is also shown. The accuracy of the present solution cannot be 
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Fic. 7. Deflection of a plate with clamped deflected edges. No load. AAw=0, 
w given at edge, w, =0 at edge. Recorded: AAw-w. Solution time: 12 hours. 1st solu- 
tion: 6 hours. Checking and eliminating mistakes: 6 hours. 


computed, since the exact solution of the differential equation is not known for these 
boundary conditions. However, one’s confidence in its precision is increased by ob- 
serving that the finite difference solution of a square plate under uniform load is in 
error less than 1% in maximum deflection when only 9 interior points are used [11], 
rather than the present 22. 


C. The equation of natural modes of a membrane, 
Aw + Aw = 0. (19) 


In problems of the vibration of two-dimensional systems, the information sought 
concerns the natural frequencies and the characteristic functions. In the equation, 
it is required to find the permissible values of \, which is the square of the frequency 
times certain physical constants. The relaxation method can be applied to this prob- 
lem in several ways. For example, a value of \ could be estimated by Rayleigh’s 
principle from an assumed deflection w, 


i al P » (w*. + wy) ; (20) 


le 
where the summation extends over all the net points. This value could then be in- 
serted into Eq. (19) in the finite difference form, 


Wi + We + ws + we — (4 — Ab?) wW = OQ, (21) 
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and the relaxation process used to reduce the Q. Periodically during the solution A 
would be reevaluated and reinserted into Eq. (21). Thus new residuals would appear, 
which could again be reduced, the process being repeated until Eqs. (20) and (21) are 
satisfied with sufficient accuracy. This method of solution, as used by Southwell [12], 
imagines that at each net point there is a force Q applied and that these forces are to 
be removed by changes of the amplitudes w, and the “frequency” X. 

Another method which seems to be superior in some respects has been worked out 
by Dr. A. Vazsonyi, and will be published shortly.* Equation (19) is written in the form 


Aw 
N= --—: (22) 


wW 


At each net point a value of the amplitude w is assumed and the corresponding values 
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Fic. 8. Lowest frequency and natural mode of quadrangular membrane. 
Second approximation. 


r 
Av+—v=0, v = 0 on edges, 
* ZvAv a 


6 = 1, a = 24, An = = .8984, A= — An. 
20? & 





First characteristic value: \ =21.562. 


( v 





Recorded: . 
- 1000 — (= 1000A*) f 








* J. Appl. Physics, 15, 598-606 (1944). 
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of Aw and A recorded. These assumed amplitudes constitute the solution to the prob- 
lem of the vibration of a membrane of non uniform mass distribution. For a uniform 
mass distribution the value of \ should be the same at each point. The values of w 
are changed to equalize the values of \. Whenever a w is changed the changes of Aw 
are immediately computed by the Laplace relaxation pattern, Fig. 3. New values of \ 
need not be computed every time. At any stage of the solution the correct X lies be- 
tween the highest and iowest value at the net points. Finally a value of \ is computed 
as an average of those at all the net points. The best theoretical way to compute 
this “average” is by the use of Eq. (20) in the form 


) wiw 


he—S (23) 


Dw 


Fig. 8 shows the solution for the lowest frequency of a quadrangular membrane. 
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Fic. 9. Second frequency and natural mode of quadrangular membrane. 
Second approximation. 


r 
Av +—v7=0, v = 0 on edges 
a 
* SvAv a « 


6 = 1, a= 24, An = = = 1.6716, = zo 





Second characteristic value \ = 40.118. 


Recorded: | v 
A ; 
|- 100 — (= 100n*) 
\ v 
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To compute higher modes, the procedure is exactly the same as outlined above 
for the lowest mode. The assumed amplitudes should of course include one or more 
nodal lines, depending upon the mode sought. All of one’s information about the vi- 
bration of membranes should be used in selecting the amplitudes, and generally it is 
convenient to make use of the orthogonal properties of the characteristic functions. 
Thus, it is always possible to use any arbitrary amplitudes w and remove the portion 
of these arising from the first mode w (obtained at the end of the first mode computa- 
tion), to obtain the amplitudes w, of the higher modes only; 
Wn = W— AU), (24) 
where 
> WW 
. aa (25) 


¢=-—=; 
> Wj 


The second mode of the quadrangular membrane is shown in Fig. 9, while Table III 
compares the results. We note the good accuracy obtained with a few points in the 
square. By the use of symmetry, the 9 point approximation for the square required 





only 10 minutes. 
TABLE III 


Characteristic values of membranes 





Square membrane; Membrane accord. to Figs. 8, 9; 























Characterization of problem lowest mode lowest mode second mode 
Order of approximation | first second | exact | first | second | first second 
Degrees of freedom of approx. 

system 9 49 — 7 15 7 15 
| | eat er | - , 
Characteristic value | 18.75 | 19.508 27? = 20.61 | 21.562 | 37.10 | 40.118 
19.739 
Number of modifications used 5 35 —— 12 35 | 16 65 
, Accuracy 5% 1.2 | 














It should be noted that the solution of a forced vibration problem is best carried 
out by the first method outlined above. 


D. Equations of the type, 
Prez t f(%, Vs Oy Par Pu) Guy = B(X, Vr Oy Oxy Pu)- 


(26) 


Equations of this type are of frequent occurrence in engineering problems, but 
because of their non-linear aspects only limited assistance is offered by conventional 
mathematical methods. For so complicated an equation, there are many possible ways 
of applying the relaxation process. Only one will be mentioned here. We obtain a 
finite difference equation by transforming the second derivatives only; 


gi + g2 + foys + foys — 2(1 + fo)eo — god = Q, 


for which the relaxation pattern is shown in Fig. 10. Thus the pattern is different 


(27) 
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for each net point and varies during the course of solution. The luxury of an investiga- 
tion of the classification of equations of the type of (26) (as to whether they are ellip- 
tic, parabolic, or hyperbolic, if such a classification is possible), of the nature of 
solutions, of the permissible boundary conditions, etc., is denied during war time by 
the urgency to get numerical results. Only casual observations have been made to 
date, and will not be discussed. However, it is certain that for f>0 [ Eq. (26) of elliptic 
type] the process of solution described is quite easy and quick to carry out. 

The domain of the problem to be solved is drawn and a 
solution is guessed at a net of points. From this set of values 
of ¢, fo and go and then the Q are computed at each point. 
By the relaxation process carried out exactly as described 
for Laplace’s equation, except that the influence coefficients 
of Fig. 10 are used instead of (—4, 1, 1, 1, 1), the residuals Fic. 10. 

Q are reduced somewhat. Before bothering to eliminate 

the Q completely, we compute new values of fo and go and hence corrected values of 
the residuals attached to each net point. This process of reduction and correction is 
continued until sufficient accuracy has been attained. Just as in previous problems, 
a finer net can be added for greater accuracy. 

As an illustration, let us consider the distribution of electric potential in the space 
between two parallel planes, one of which has, standing normal to it, a right circular, 
cylindrical post with a hemispherical top. Fig. 11 shows the problem and the solution. 
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Fic. 11. Axially symmetric electric potential distribution. 
yp 1 &» e 
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If cylindrical coordinates are used with the z-axis along the axis of the post and the 
origin at the base of the post, the potential is described by the equation 


0? 1 @ re) 
S+— = (r2)~0 (28) 
02? r Or or 





ec 


with the boundary conditions ¢=0 on the base plane and post, and ¢=100 on the 
parallel plane. Again, several alternative procedures are possible. For example, Eq. 


(28) can be written 


1 
zz + Cnr = —_——" &, (29) 
r 
and the solution can be carried out as described above for the general case. 
In the case of the present problem it was decided to make the substitution 
y = logr, (30) 
thus converting Eq. (28) into 
1 
Pz: + =a 0. (31) 
y? 
Hence the finite difference equation becomes 
¢1 + os + Bigs + gs) — 2(1 + B)go = Q, (32) 


where B= 62?/r?5y?, dy and 6z being the net spacing in the y and 2z directions, re- 
spectively. The net spacing was chosen so that 6z=4éy. In this way sufficient points 
appear where they are needed, i.e., near the post. In Fig. 12 the transformed plane is 
shown, together with the influence coefficients at the top of each column of points. 


E. Equations of the type, 
) o 
— (ugz) + — (uey) = B(%, ¥, 9, Px Oy), (33) 
Ox oy 


where p=p (x, ¥, , Gz, Gy). This equation is of very general occurrence in physical 
problems. For example, in the case of a soap film with large deflections w under con- 


stant excess pressure p, we have [14] 


= constant; (34) 





g= w, t= 


1 p 

’ su 
{1+ wi + wi}? 2 
in the case of the plane irrotational flow of a compressible fluid with velocity poten- 


tial g, we have 





. y-1 0%) 1 
w= p = density = po41 — “. , g = 0; (35) 
1 a)yy-—1 


in the case of a steady magnetic field (g = magnetic potential) in a non current carry- 
ing medium, we have 
= . “qe = 2 2 
= magnetic permeability = f(¢; + ¢)), (36) 


where the function / is given by experimental data on the material. 
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The flow of oil in the film lubrication of a bearing, when the viscosity is assumed 
independent of temperature, is described by an equation of this type [see Eq. (39)]. 
In general, physical problems which lead to Laplace’s equation when physical param- 
eters are constant, give equations similar to (33) when the physical properties vary. 
The properties may be known functions of x, y because the material is non-homo- 
geneous, or of ¢, @z, ¢y, etc., because of the nature of the material itself. 
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Fic. 12. Transformed y, s-plane of axially symmetric electric potential problem. Recorded: ¢. 


For the numerical solution of these problems it has been found convenient to 


use the form 
Grz + yg = Em (log h) zPz — (log BH) yy (37) 
In finite differences this becomes 


01 + ¢2 + 3 + Gs — 40 
= god? + 3{(¢1 — ¢s)(log ws — log wr) + (ga — 2)(log wz — log ws)} +0. (38) 


The relaxation pattern used is that used for Laplace’s equation. The Q is to be elimi- 
nated. The variable terms on the right of Eq. (38) are computed periodically as cor- 
rections. This method is well adapted to the equation as long as Eq. (33) remains of 
elliptic type. The “correction” terms on the right contain ¢2 and gy through the 
relation (34). When these variations are such that Eq. (33) becomes hyperbolic, the 
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relaxation process becomes confusing unless a great deal of physical knowledge is 
available to assist the computer. 

The solution of a problem of this kind is shown on Fig. 13, which shows the shape 
of a soap film with large deflections (Eq. 34). The maximum deflection is 100 units, 
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Fic. 13. Soap film with large deflections. 
a a 1 
=~ (uwz) ts (uwy) = 0, where p= a+ wo +n 
Recorded: deflection w. Time to correct small deflection solutioa (shown by broken 


lines): 3 hrs. 


compared with the largest dimension of 160 units for the xy projection of the bound- 
ary. We note the deflection as given by Laplace’s equation and show it for comparison 
by dotted contours. In solving this problem, the Laplace equation solution was taken 
as the first approximation. Correction terms had to be computed only twice. This 


solution required 3 hours. 


F. A more general type of equation which arises in the flow of oil in a bearing has 
been solved by Christopherson [13]. He solves the system of equations 


oF? ef ofH* oP OH 
pee ah), ape ory oe, ro 
0ELM oé Oni. M On 0g 
oT H* oF oT HH? oP M oP 
2a Pie a a Pea a 
dé M do On M on H? 0g 


where H is a given function of £, 7 and M is a given function of T, P; P is the pressure 
in the lubricant, T is the temperature of the lubricant. 
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G. Another very general system of non linear, integral, differential, difference equa- 
tions arises in the problem of the thermal equilibrium of a nest of N coaxial, conduct- 
ing cylinders of length 2b, between which a hot gas flows. The temperature T of 
the cylinders varies with the number of the cylinder and with the axial position z 
of the point under consideration. The temperature of the gas varies with the number 
of the stream, and the position x. The equations to be solved are 


Silx, n)T*(x, n) + of {T4(y, nm + 1)fo(x, m, y) + Ty, n — 1)fs(x, n, y) }dy 
—b 


0°T (x, n) 
+ b{T.(x, n) + T.(x, n — 1) — 2T(x, n)} + i . =0, (41) 
x 


2n+ 1 


n 


OT, (x, n) 


Ox 








n+1 
4 {7 (2, n) + —— T(x, + 1) —- T (2, mba = 0, (42) 
n 


where a, b, c, d are constants 

Si, fo, fg are known functions of the variables indicated. The solution of these equa- 
tions, which required about 30 hours, is shown in Fig. 14. The temperature contours 
show only the distribution of cylinder temperatures. This analysis was actually used 
as the basis for redesign of the instrument involved. 
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Fic. 14. Thermal equilibrium of a nest of cylinders in a hot gas stream including 
radiation, conduction and convection. Gas temperature: 1700° F abs. Surrounding 
duct temperature: 1500° F abs. Section of nest of cylinders showing cylinder temper- 


ature distribution. 





4. General remarks on finite difference approximations. In this paper only the 
simplest approximations for rectangular nets of points have been mentioned. There 
is no reason why higher order approximations cannot be used, or why triangular or 
other point arrangements should not be considered. Indeed both of these possibilities 
have been used [9, 10, 13]. There is no sure way of deciding at the present time just 
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what net should be used with a given problem. The ultimate requirement is a suff- 
ciently accurate answer in the shortest possible time. The author’s experience to 
date favors the simplest possible formula, the simplest possible relaxation pattern, 
with accuracy obtained by extra points on a finer net where needed. Southwell and 
Christopherson [9, 10, 13] have occasionally found other formulas to be advantageous. 
Until a person is thoroughly acquainted with the details of the several processes 
of solution, the “bookkeeping,” it is impossible to judge the relative merits in the 
matter of time required to reach a given accuracy. 

This is particularly true in the matter of boundary conditions. Every new com- 
puter “discovers” the possibility of deriving special formulas to apply to points near 
(or on) the boundary of the domain, especially when the boundary itself wanders 
among the net points. In many cases such formulas are admirable from the point of 
view of accuracy for a given net size, but fail miserably when compared with linear 
interpolation (or extrapolation) together with a somewhat finer net (used locally). 
By far the best general procedure for fitting “queer” boundaries or boundary condi- 
tions is to sketch a graph for the last few points approaching the boundary and thus 
compute the boundary point graphically. 

5. Conclusions. This paper stresses particularly the practical aspects. The finite 
difference approximations to partial differential equations are well known; only the 
detailed steps in carrying out the solutions are not yet general knowledge, and these 
are discussed step by step. The relaxation method, first conceived by R. V. Southwell, 
is the underlying process in the solution of most of the problems discussed. Once the 
basic idea is grasped by actually solving a problem, it is capable of enormous extension 
to a great many kinds of problems with only the most meagre knowledge of the 
current methods (if any) for solving them. From an engineering point of view this is 
of enormous importance, because the practical man would like to do something better 
than guess, and yet he cannot afford the time required to become versed in analytical 
procedures, procedures which too often cannot supply a numerical answer to the real 
physical problem with reasonable accuracy and speed. 
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STUDIES IN OPTICS 


I. General Coordinates for Optical Systems 
with Central or Axial Symmetry* 


BY 
M. HERZBERGER 
Communication No. 960 from the Kodak Research Laboratories 


In previous papers!” the author has proposed an approach to gedémetrical optics 
different from that developed by Hamilton and his successors. 

The purpose of the present paper is to generalize the formulas in these papers, 
and to find the most general treatment of systems with central (point) symmetry and 
with axial symmetry. By leaving the coordinates general, subject only to the sym- 
metry conditions of the problem, we retain the symmetry in the formulas up to the 
point where we desire to draw conclusions for a special problem. We can then intro- 
duce special coordinates adapted to the problem in question, and find the particuiar 
answers. 

The fundamental-invariants of geometrical optics show no preference for either 
object or image side, nor for point or angle coordinates as variables. The different 
approaches suggested by Hamilton, as well as the direct approach just mentioned, are 
special cases of the methods developed here corresponding to special choices of co- 
ordinates. Several different choices of coordinates will be given as examples. 

The fundamental formulas (A, B, B’, C below) are based only on symmetry con- 
ditions and on the validity of the Lagrange invariant (A). They are therefore not 
restricted to optical problems,* but are also valid for problems in mechanics, hydro- 
dynamics, and electron optics. 

1. Ray tracing formulas, the Lagrangian invariant. Let us assume a ray travers- 
ing a number of optical media with refractive indexes m, m2, m3, - - - , nm’. Let a(x, y, 2), 
a’(x’, y’, 2’) be a vector from an arbitrary origin to a point on the object and image 
rays, respectively. Let ax(xz, ys, 2) be the vector from the same origin to the inter- 
section point of the ray with the &th surface. Let S¢,241(£:.241, Me,e+1) [k,k41) be a vector 
along the ray in the medium between kth and (k+1)th surface, a vector of length 
equal to the refractive index m;,441 of the medium. 

Let 0, be a vector perpendicular to the kth surface at the intersection point. Its 
length may remain arbitrary, for the moment. The refraction law then reads: 


Sk,k+1 X O, = Sr-1,k X Ox, (1) 


where the multiplication sign indicates vector multiplication. Equation (1) shows 
that S:.441—S-1,4 has the direction of the surface normal o,, or 





* Received Jan. 24, 1944. 
1M. Herzberger, Direct methods in geometrical optics, Trans. Am. Math. Soc., 53, 218-229 (1943). 


2 M. Herzberger, A direct image error theory, Quarterly of Applied Mathematics, 1, 69-77 (1943). 

3 For the connection of the Lagrangian invariant with different branches of mathematics and physics, 
see M. Hetzberger, Theory of transversal curves and the connections between the calculus of variations and the 
theory of partial differential equations, Proc. Nat. Acad. Sci. U.S.A., 24, 466-473 (1938). 
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Sk,k+1 — Sx—-1,4 = PxOx. (2) 


We now can describe the path of the ray through the system by means of the vector 


equations 
a + As, Sic = S + $10), @. = a; + A2Si2, 


. * (3) 


ai 


s’ -_ S,-1,» + ,0,, a’ = a, a d's’. 


The geometrical significance of \ and @ can be seen by multiplying (3) scalarly by 
S441 and 0,, respectively, keeping in mind that, by definition, 


Stati = Meat (3a) 

We then find that 
Nene = (Ae — Acs) *Se—1,4/Me—-1,k; (3b) 
bk = (Sk,n41 — Sk1,x)*Ox/0}, (3c) 


i.e., the \’s are proportional to the distance between the two surfaces along the ray, 
and the ¢’s are proportional to what might be called the power of the surface for the 
individual ray. 

Since Eqs. (3) are valid for every ray, we now consider a two-dimensional mani- 
fold of rays, i.e., we assume the a’s, 0’s, and s’s to be vector functions of two variables 
t; and ft. From the definition of s,,4,1 and 0%, we find 


Si,e¢1° (OSz,n41/0t,) = 0, 0,.- (da,/dt,) = 0, (u = 1, 2). (4) 


We now differentiate (3) with respect to 4; and multiply scalarly by Os:,441/0/ and 
0a,/Ot, respectively. Then we differentiate with respect to f and multiply scalarly 
by Os¢,441/0h, and 0a;/0t,, respectively. Subtraction of the two sets of equations yields 
the “Lagrangian invariant”: 




















> =) = a a a 

| Ot, oat = Ot, Oty a te | Ot, at; | (A) 
da OS | | 0@; OS1,2 | da’ as’ 
at, oh | | Oh dh | | at, ote 


This formula was introduced by Lagrange in his astronomical investigations. It is 
known by the name of the Lagrangian bracket in the theory of partial differential 
equations. Herzberger* used it in his theory of transversal curves, as the starting 
point. Let us now see what conclusions can be drawn if the system in question fulfills 
certain conditions of symmetry. 

2. Centrally symmetric systems. In this case all refracting surfaces are concentric 
spheres with radii rm, ---, f,. It is therefore appropriate to consider the common 
center as the coordinate origin, and to choose concentric spheres as the object sur- 
face a and the image surface a’. All the surface normals pass through the common 
center. We shall give them the length r; from center to surface, with a positive sign 
if the surface is convex towards the incident light. In other words, we make 0; =a. 
Under these conditions, Eqs. (3) become 
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a, =a-+ds, Si2 = S + $181, 
ew (5) 
s’ = Sn—1,n + PnAn, a’ = a, + ds’. 


From Eqs. (3) we can find an invariant vector, namely, 
aXs=a:Xs =a, X82 =-:: =a’ XS =p. (6) 
Therefore, in a concentric system, both object and image rays lie in a plane through 
the center, and the optical length p of the perpendicular from the center (the length 
of the invariant vector) remains constant. 
Equations (4) can be combined into 


a’ = aa + bs, s’=ca+ds, where ad — bc = 1. (B) 


The invariant relation, ad —bc =1, is found by substituting (5) in (6). It is possible to 
calculate a, 6, c, and d as functions of the \’s and @’s, if we use “Gaussian brackets.”* 


We find 


a 


[¢1, A12, en d’], b = [n, 91, ee Te r’], (7) 
/ 
[¢1, 12, nea fk onl, d sated [n, $1, "oy oe on |. 


In the case of central symmetry, the ¢’s and X’s, and therefore a, b, c, and d, can 
be considered as functions of a single variable, and p can be taken as this variable, 
as shown in (9) and (10). Now 

(@x°Sz,n41)? + (Qe X Skn41)? = St n+1) (8) 


or 
O4-Skis1 = Vries, — p?. (8a) 


Il 


C 


Thus we find from (3) that 
1 


< a aac i. ae 

ne ee 2 ao see on 

Ate = a [s Tr 1Me e+ p V Meet p | 
ky, k+1 


1 2 2 
Oe ee ee) 
TKe+INk,k+1 TENK k+1 


Nk jk+1 








and that 


1 33 
n= 5 [Vln — F ~ Volt ~ Pl (10) 
k 


Equations (B) correspond to the direct equations of our theory. A more general 
representation is given by choosing two arbitrary vector functions, 1 and m, in terms 
of which the object and image vectors can be expressed. Equations (B) may then be 


written, 
a = al + dm, a’ = a/1+ b/m, 


s = cl + dm, s’ = cl + d’m; (B’) 
ad — bc = a'd' — b'c’. 


4 Invented by L. Euler in 1776. See M. Herzberger, Gaussian optics and Gaussian brackets, J. Opt. 
Soc. Amer. 33, 651-655 (1943). 
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a,b, c,d, anda’, b’,c’, d’ are still functions of a single variable, which can be taken as 
x = (1 X m)? = I’m? — (1-m)*. (11) 


It is obvious that under these conditions the image formation described by (11) fulfills 
all the conditions mentioned, and that the last condition in (B) is equivalent to the 
validity of the invariant (6). 

Formulas (B) are a special case of Eqs. (B’) if we choose 1=a, m=s, which corre- 


sponds to a=d=1, b=c=0. 
The connection between a, b, c, d and a’, b’, c’, d’ becomes clearer if we introduce 


some auxiliary angles. Let us write 
<(a,m)=a, <X(Lay=8, <“C(m=y  <£1,s) = 4, (12) 
< (1, m) = y, X (a, s) =a, X (a’, s’) =’. 
From (B), if we write r for the absolute value of vector a, we find that 
a = (rm/zx) sin a, b = (rl/x) sin B, 


c = (nm/zx) sin y, d = (nl/zx) sin 6, 


(13) 


with analogous expressions for the primed quantities. Moreover, we find from (12) 


that 
atBp=ytb=P=ao't+pP=7 +8, 


a-y=6-fB=a, a —-y =i -p =, (14) 
where o and a’, according to (9), are connected by 
p = nrsino = n’r' sino’, x = Im sin ¢. (15) 
If we write 
ad — bc = ad’ — b'c’ = D, (16) 


we have finally 
p = Dr. 

Thus, for a given system of coordinates 1, m, and ¢, and given object and image 
spheres (radii r and r’), we can calculate all the functions a, b, c, d, a’, b’, c’, d’, if 
only one of them is given on each side. For instance, let us assume y and yy’ to be 
given. We then find o and a’ from (15), and obtain 


a= : =y-—o-y, =y, b=y-Yy 
o+y¥ B=y-c-y Y=y7 Y¥-Y (17) 


U 


d=ad+y7, B=v-d-y, =r, 8 =v-7. 


II 


Thus, according to (13), we determine a, 5, c, d, a’, b’, c’, and d’. 
Let us now consider some special cases. 
a) The direct method. We choose 1=a, m=s. This means that 


y=ao, B=y7=0, am=z=j=oe, 
a=d=1, b=c=0, (18) 


aad’ t+y, P=a-d—-y, Yu, Feo-y, 


where mr sin g=n’r' sing’. 
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b) Hamilton's point coordinates. We choose 1=a, m=a’, or a=1, b=0, a’=0, 
b’ =1. This means that 


rn=y, = 0, =y-— a, 6= 0, 
a=y B y=¥ (16) 
a’ = 0, p’ = yf, y' = — oa’, f=yt+o’, 
or, since l=7, m=r’, 
a= 1, b = 0, c = [nr' sin (¥ — o)|/x, d= p/x, 17’) 
7 
a’ = 0, b’ = 1, c! = p/x, d’ = [n’r sin (W + o’)|/z, 
where 
mr sino = n’r’ sin o’ = fp, rr’ siny = =. (18’) 
We note especially that c’—d =0. 
c) Hamilton's angle characteristic. We choose 1=s, m=s’, or c=1, d=0, c’=0, 
d=1. We find from (B’) that 
ad — bc = a’d’ — b'c’ (19) 
or —b=a’'=D=p)/rm. For the auxiliary angles, we obtain 
r=ot+y, =-«, oe ae Tn 
c v B ,=¥ (20) 
a’ =o’, B=y-’, y’ = 0, i =y, 
or 
a =(n'r/x) sin (YW +0), b= —a’ = — p/n, Bb = (n?/x) sin (¥ — 0’), (21) 
where 


x = nn’ sin y, p = arsino = n’r’ sin o’. 


d) We take as coordinates the intersection points of the ray with two spheres in 
the object space: the object sphere (radius r), and a second sphere (radius R, and A, 
the vector to the intersection point on this sphere). We define 


1 =a, m= A, 
and 
a= a, a’ = a’'a+ 0’A, (22) 
s=-ca+dA, s'=ca+d’A, 
where a’d’—b’c’ =c and 
c = [nR sin (¥ — o)]/e = — [nr sin o]/x = — p/n. (23) 
We find that 
sin o’ = [nr sin o]/n’r’, sin (Wy — oc) = [n sin o]/nR, 
and, finally that 
a’ = [r’R sin (o’ + y’)]/z, c’ = [n’R sin y']/z, (24) 


b’ = [rr’ sin (W — o’ — y’)]/x, d' = [n’r sin (Wy — 7’)]/x. 


These are the equations for the image formation. 
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3. Systems with rotational symmetry. Let us now assume that the system has sym- 
metry only around an axis, the unit vector along which we shall designate by k. 
In this case, all the surface normals intersect the axis, and we shall give to the vector 
0; in (3) the length r;,, which is the distance along the normal between its inter- 
section points with the axis and the surface. 

We now project all the vectors on a plane perpendicular to the axis, and define 
the projected vectors b,; and tx,x41 by the equations 


a, = b, + 2k, Sicgr = teega Ht Seegik, oO. = bi + (2: — 2we)K, (25) 


where zy is the quantity known in geometry as the subnormal, and 





2 2 2 zt 
tke = V ness — ets — M+ 


x,441 is the (optical) cosine of the angle between the ray and the axis. 

Let us now assume that object and image origins lie on two planes perpendicular 
to the axis. We can then replace all the vectors in (3) and (A) by their projections 
in these planes, and find, instead of (A), for a two-dimensional manifold of rays 
(parameters 4, /2), 











db sat | | db’ at’ 
ot ot ot Ot, | 
| in ‘ ‘ a at’ j — 
| At, at, | | dtr ate | 
and instead of (B), 
b’ = a’b + U't, t' = b+ dt, (27) 


where b’ Xt’=b Xt and therefore a’d’ —b’c’ =1. The functions a’, 6’, c’, d’ are given 
by formula (7), where ¢@ and \ have the same meaning as before. 

Moreover, a’, 6’, c’, d’ are no longer functions of a single variable, but are func- 
tions of the three symmetric functions and b and t, namely, 


u= b?, v= bD-+t, w= t?. (28) 


Equation (27) corresponds to the formulas of the direct image error theory. The 
most general choice of coordinates might be described as follows (t and m, as well 
as the other vectors, lie in a plane perpendicular to the axis): let 


b = al + dm, b’ = a‘1+ b’m, 


(29) 
t= c+ dm, t’ = cl + d'm, 
where 
bxXt=b’' xt’, 
or 
ad — bc = a’d’ — Bc’. (30) 


Let us assume that a, b, c, and d are functions of the symmetric functions of 


land m; that is, 
@ = |, v=1-m, w= m?, (31) 


b, t, b’, t’ must fulfill Eq. (26), if we set 4 and # alternatively equal to u, v, and w. 
Thus, we find the following equations for a, b,c, and d. Let us write 
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au + bv, C = cu + dp, 
av + bw, D= c+ dw, 


A 
B 


(32) 


H 


and introduce the abbreviation A for the difference of an expression in the object and 
image spaces. An easy computation then gives 








| 04 94 | | OB OB da da | | db ab] 
Ot; ate | Ot, ate | Ot; te | dt, te 
A + + : a = 0, (C) 
Oc Oc | Od odd | eC @ | | dD aD 
Ot; Oto Ot; Oto | Ot; Ote | | Ot; Ote 


A(ad — bc) = 0. 
These are the necessary and sufficient conditions that (29) describe an axially sym- 
metric image formation. We repeat that Eqs. (29) and (C) demand only the validity 
of Lagrange’s invariant (26), and axial symmetry. Their application is therefore not 


restricted to optical problems. 
Let us now again investigate what forms the fundamental formulas assume for 


special choices of coordinates. 
a) Hamilton's (Bruns’) point characteristic. Hamilton (Bruns) chose as variables 


the coordinates of a point in the object and image spaces. That corresponds to taking 
l=b, m=b’. Equations (29) become 
b = b, b’ = b’, 














(33) 
t = cb + db’, t' = cb+d'b’, 
or 
a =D’ = 1, b=a =0. 
These are conditions for the coefficients. Equations (32) now become 
A = 4, C = cu + dz, A’ = 2, C’ = cut d’2, (34) 
B=», D= c+ dw, B’ = w, D=cvr1+d'w, 
and we find instead of (C) that 
d=-d¢, (35a) 
| du Ou | | ov av | 4 | dd’ ad’ 
dh a | | ah Of | | dn a | | dh abs : 
1+ | | + | = 0. (35b) 
0c Oc | Od dd | Ov Ov Ow Ow | 
Sas, coca ee eeeeet min fate 
Ot; Ole Ot; Ots Ot; Ot» Ot; Ote 


Equation (35b) stands for three equations, which we can obtain by replacing 4; 
and #, in (35) by u and v, u and w, v and w, respectively. This yields 


0c = 0d Oc’ Oc dd’ od Oc’ od’ 
=——-», — = —_ —- —_ - (35c) 


———ia— Sp — = 


Ov Ou Ou Ow Ou ow ow Ov 
Equations (35a) and (C), when integrated, lead to a function V(u, v, w) such that 


OV OV OV { 
Poe d aie. uniadiains 


re] 
c= 2—» c= -—») = , d’ = — 2—-- (36) 
Ou Ov Ov Ow 
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V is the characteristic function of Hamilton, the “eiconal” of Bruns. Formulas 
(36) agree with Hamilton’s formulas, except that he used 1*/2 and m?/2 as variables. 
Our choice of variables simplifies the form of the general formulas (C). 

b) The angle characteristic. Hamilton chose as coordinates the direction cosines 
of the object and image rays. This means that 1=t, m=t, or 


b = at + dt, b’ = at + b't’, 
(37) 
t=t, =f’, 
or 
c= 1, d = 0, c’ = 0, d’ = 1, (38) 
Equations (29) now give 
C=u, C’ = 9, D=», D’ = w, (39) 
and (C) becomes 
b+ a’ =0, 
0b da da’ da ab’ 0b da’ —s 0b’ 
——-— = —-, —-— = — —-—= —-—+— (40) 
Ou dav Ou Ow Ou Ow Ow dv 


Equation (4) is solved if we introduce the angle characteristic T(u, v, w), and set 


1 aT aT aT 1 oT 
‘=—-—, be— Y¥u-——- (41) 


oe SS ses e¢= 


2 Ou Ov Ov 2 dw 


We see that this also agrees with Hamilton's theory. 
c) The direct method. In the papers mentioned,':? we took as variables the object 
point and the direction of the object ray, i.e., we chose |= b and m=t. This gives 


b =b, b’ = a’b + Ut, 





(42) 
t=t, t' = b+ ¢d't. 
That is, we put 
a=d= il, b=c=0. (43) 
Equations (27) then give 
A=4, C= 9, B=z2, D= »w, (44) 
and Eggs. (C) give 
a'd’ — We’ = 1, 
| 0A’ @A’ OB’ AaB’ | da’ Aa’ | | 0b’ = ab’ 
| Ot, be Ot, Ate Ot, Ale | Ot, te 
+ + + = 0. (45 
dc’ ac’ ad’ ad’ | ac’ ac’ | aD’ aD’ | 
At, dt, Ot, te | | dt; Ot, | | dt; ate | 


If we denote the sum of the four determinants in (45) by J’ when 4, =4, 4 =v, by 
II’ when 4;=4, 4 =w, and by III’ when t;=0, t2=w, we can write (45) in the form 


a’d’ — bc’ = 1, =I] =III' =0. (46) 
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Again, if we disregard the fact that the variables u and w differ by a factor of 
two from those used previously'? we find Eqs. (46) to be identical with those given 


before. 


d) Object and stop coordinates. To analyze image errors, we investigate the manner 
in which the image ray changes with the position of the object point and the position 
of the intersection with the diaphragm plane, for which we frequently substitute the 
entrance pupil of the system. If we choose these as the coordinates of the ray, assum- 
ing that the distance between object and entrance pupils is equal to k, we find that 


b =b, b’ = a'b + b’b,, 
t=7(b—b,), t=cb+db,, 


where 
n Nn 


JVE+b—b) VEtu-btw 











From (48) we conclude that 
1 1 n 
4 ee ee ae Te 2 eee a + 
- 2 V(k? + u — 20+ w)' 


Thus Eg. (29) gives 


A=4, B=z1, C = y(u — v), D=y(v — w). 


The fundamental equations (C) become 
n 


a’d’ —p’ ‘= Te 
VR+u—%w+w 





is 


n 
’=-]!' =]1I' =— --- . 
2 (k? + u — 20+ w)?!? 





(47) 


(48) 


(49) 


(50) 


(S1) 
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LATERAL BENDING OF SYMMETRICALLY 
LOADED CONICAL DISCS* 


BY 


K. E. BISSHOPP 
Fairbanks, Morse & Company, Beloit, Wisconsin** 


1. Introduction. The general theory of lateral bending for thin circular plates of 
variable thickness is given in Timoshenko’s book,’ “Theory of plates and shells,” 
where also may be found numerous references to the literature of the subject. Of 
particular interest is a reference to Féppl who indicated the analogy existing between 
the rotating disc problem and that of lateral bending in a circular plate of variable 
thickness. Comparison of the solution for the rotating conical disc problem with the 
corresponding one for lateral bending shows that the basic differential equations in- 
volved, and the expressions for the stresses, are analogous. Therefore, previously de- 
scribed methods? for obtaining solutions of the former problem in terms of hyper- 
geometric functions are applicable to the latter problem. It will appear later that the 
special type of hypergeometric differential equation associated with the lateral bend- 
ing problem has solutions which give the stress coefficients with less labor than in the 
case of the rotating conical disc. 

The stress coefficients have been arranged conveniently for numerical calculation 
of conical discs, which are component parts of a wide variety of engineering struc- 
tures. The head of a large poppet valve provides a particular example where the 
principal stress member can be approximated by a system of incomplete conical 
discs. In order to illustrate an application of the theory, stress coefficients for conical 
discs subject to lateral bending as well as for rotating conical discs will be used to 
estimate stress distributions in a steel valve head of constant weight and various 
proportions. Since the coefficients are obtained from solutions of differential equa- 
tions for thin discs, the approximate method breaks down in the neighborhood of the 
valve stem. These limitations have little effect near the periphery, which makes it 
possible to calculate valve proportions corresponding to approximately uniform stress 
distribution throughout the head. The description of the illustrative example at the 
end of the paper explains the method of calculation in detail. 

2. Derivation of differential equation. Let M, and M, denote radial and tangential 
bending moments per unit length acting on an element of a circular plate at distance 
r from the center; then if Q is the corresponding circumferential shearing force per 
unit length, the equation of equilibrium is 


M, + rdM,/dr — M, = — Qr. (1) 


If w denotes downward deflection of the middle surface, then 


* Received Jan. 29, 1944. 

** Now at Armour Research Foundation, Chicago, III. 

1 Timoshenko, Theory of plates and shells, McGraw-Hill, 1st Edition 1940, Art. 54, p. 282. 

2 K. E. Bisshopp, Stress coefficients for rotating discs of conical profile, Journal of Applied Mechanics, 
Vol. 11, No. 1, March 1944, pp. A1-A9. 
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d?w o o dw dy o 
— D(— = D(*+— 9), 
dr? r a. dr r 


1 ~l oe w 
M, = at —~ => 
a ae 


where ¢@ is Poisson’s ratio, p= —dw/dr and D=Eh*/12(1—07), E being Young’s 
modulus and & the thickness of the plate; D is called the flexural rigidity. 
In the case of a conical profile, /# is a limear function of 7, so that substitution 


from Eqs. (2) in Eq. (1) gives 
d {dg 7) dD {de ¢g 
p—(4+*)+—(F+0*)- —Q. (3) 


dr \dr r dr \dr r 


M, 


ll 


(2) 


In order to reduce this equation to non-dimensional form, we introduce the radius R 
the knife edge of the disc and the thickness hy at the center. If r/R=x, then 
=h)(1—x), and Eq. (3) becomes . 


d%y 1 3 \de /1 30 120R*%(1 — 0?) | 
—+(—-—_)" -(4+=)e-- - (4) 
dx? x 1— x/ dx x? «(1 — x) Eh3(1 — x)* 








For any particular type of symmetrical loading the shearing force Q is a function 
of x alone. The maximum radial and tangential bending stresses S, and S; are ob- 
tained from the general solution of Eq. (4) with the aid of Eqs. (2) and the relations 


S, = 6M,/h’, S; = 6M,/h’. (4a) 


The problem of a conical disc supporting a concentrated vertical load P at the 
center has some interesting practical applications. In thiscaseQ=P/2rr=P/27Rx, 
and Eq. (4) becomes 


—-+ 3 = — —— (5 
va wEhi(1 — x)? 





d*y dy 1-—x 6PR(1 — 0%) | 
ax(1 — x) + (1 — 4) 2 — ( an) ( 
dx? dx 


x 


It can be verified by substitution that a particular integral of Eq. (5) is 


2PR(1 + a) — 3c 1 1 
_ (= +-+—). (6) 
wEh(1 — 3c) \(1 — x)? x 1—z 








$3\ x 


The auxiliary equation, the solutions of which are independent of the type of loading, 
is obtained by setting the right hand side of Eq. (5) equal to zero. After making the 
substitution g=xF, we obtain 

a°F a 


d 
e(1 — x) —-+ 3(1 — 2x) — 3(1+ 0)F = 0, (7) 
dx? d 





x 
which is recognized to be of hypergeometric type. 
3. Complementary functions. Equation (7) is of the form 


9 
é 


d°*F dF 
“(1 — x) —-+ [c — (a +6 + 1)2] — — abF = 0, 
dx? dx 
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where c=3, a+b=5, and ab=3(1+0). The first solution can be represented by a 
power series; the integral exponent difference* 1—c = —2 shows that the second solu- 
tion contains a logarithm.* In the notation of the hypergeometric function, the first 
solution is 
ab a(a + 1)b(b + 1) 
F,(x) = F(a, b,c, 2) = 1+—x#+ ePf+eee, 8 
l-c 1-2-c(c + 1) (8) 
which converges absolutely and uniformly when | x| <1. The asymptotic behavior of 
the hypergeometric function in the neigborhood of its poles is given by* 


T(o)P(a +b —c) 
T'(a)I'(6) 


whenever c—a—b is an integer less than zero, ['(z) being the well known gamma 
function. Thus F;(x) has a second order singularity at x =1 such that, 
r(3)P(2) 2(1 — x)~* sin ax 
F(x) ~ ———-(1-2)* = 
21 I'(a)T(d) a(a — 1)(a — 2)(a — 3)(a — 4) 

which may be used to approximate the function for values of x near unity. The 

presence of singularities of lower order in the remainder term for F(x) makes this 

method unsuitable for accurate numerical work. Better approximations for similar 

functions with second and third order singularities are given in Ref. (2). 

The logarithmic solution’ of Eq. (7) ist 


(ab — 4)(ab — 6) 1 ab—6 





(1 scl x) siete (9) 





F(a, b, ¢, x) ~ 


zl 


(10) 




















F(x) = — F,(x) log, x + — — — g(x), (11) 
2 + 
where 
re > (n+a— 3)---(a— 2)(n+ 5 — 3)--- (6 — 2) ii, (12) 
na? n'(n — 2)! 
and 
&, =¥(n—3+a)+¥(n —3+ 5) — ¥(n) — ¥(n — 2) 
1 1 1 
“joa * gaged bad” Ree 
eo cee ee eee - 
2 n 2 n—2 


The principal part of expansion (11) shows that F,(x) has a second order singularity 

at the origin. The nature of the singularity at x =1 can be recognized by observing 
’ Whittaker and Watson, A course of modern analysis, Cambridge, England, 4th Edition, 1927, 

p. 198. 

* When ¢ =1/3 both solutions can be expressed in terms of rational algebraic functions. 

‘ Titchmarsh, Theory of functions, Oxford, England, 1932, p. 224. 

‘’ Forsyth, Theory of differential equations, Cambridge, England, 1902, vol. 4, part 3, p. 147. 

t The numerical value of ¢ is used since it is independent of Poisson's ratio ¢. 
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the limiting form of the mth term of g(x) which is proportional to that of Fi(x) pro- 
vided lim,... ®, remains finite. That this is the case can be shown with the aid of the 
logarithmic derivative® of the Gamma function, which gives 


1 1 1 
lim @, = 2( ; + + -) — 2y — 2¥(a) — x cot az, (14) 
a— a 





—1 a 


n— 2 


where V(a) =I'’(a+1)/I'(a+1) and y is Euler’s constant. Thus F,(x) has a second 
order singularity at x =1 of magnitude 


iti « teh =/\ — (15) 


Fe) | n— © vT 





The slow convergence of the power series near the singularities of F,(x) and F2(x) 
makes numerical evaluation of the stress coefficients for all values of x between zero 
and unity exceedingly difficult, in spite of available asymptotic approximations. One 
scheme for removing this difficulty would be to construct from the transformed differ- 
ential equation (7) two new solutions of argument 1—x and combine them linearly 
with F(x) and F,(x), as described in Ref. (2). This differential equation is invariant 
under transformation by 1—x, which brings about added convenience of calculation; 
however, considerable further reduction in computation can be accomplished by ex- 
pressing F;(x) and F:(x) in terms of symmetrical hypergeometric functions. 

4. Solutions in terms of even and odd functions. Whenever 2c=a+6+1, which 
condition is satisfied by Eq. (7), the transformation (1 —2x)*=¢ reduces the standard 
form of the hypergeometric equation to 


d°F 1 a b dF ab 
#1 — 2) +[>-(G+5+1)4] ——F=0. (16) 
dt? 2 2 ys dt 4 








The solutions of this equation as functions of x are’ 
F{}a, $b, 3, (1 — 2x)?} =Gi(x), (1 — 2x)F{4(a + 1), 3(6 +1), % (1 — 22)"} = Ga(2). 


This shows that G,(x) =Gi(1—x) and G2(x) = —G2(1—x). Since only functions of x are 


involved, 
G(x) = CF (x) + C.F 2(x), (17) 


G2(x) = D,F\(x) + DF (x), 


where Ci, C2, Di, Dz are constants; G,(x) and G2(x) are respectively even and odd 
relative to the point x=}. The series for Gi(x) and G2(x) are very convenient for 
computation when .25 <x $.50, while those for F\(x) and F:(x) are equally so when 
0 <x S.25. Since the G’s are symmetrical it is necessary to compute only one half as 
many fundamental values for constructing tables of stress coefficients as would be 
required with the F’s. From this point on therefore, the F’s are subordinated to the 
role of “helping functions,” while the G’s form the basis of all subsequent calculations. 

Returning to Eq. (17), we employ the familiar method of comparison of singulari- 
ties for evaluation of the linear factors. It is apparent from the character of the F’s 
that the G’s have second order singularities at zero and unity whose values may be 
deduced from Eq. (9). After some reduction, we obtain 


6 Ref. 3, p. 246. 
7 Ref. 3, p. 297, Example 7. 
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_ V x T'(a)T(d) E ii sin “| 
321 ($a)T' ($b) n+ © T 
Cz = Vx/[16T(3a)T(35)], 
Vx T(a)T (0) ‘ sin ar (18) 
D,= - E — lim 4, =, 
647 {4(a + 1) } 17 {3(6 + 1)} 20 ~ 
D, = Vx/[32r {3a + 1)}T{3(6 + 1)}]. 


The functions G;(x), Ge(x), and their derivatives are tabulated in Table 1. 

Since the F’s and G’s are linearly dependent, xG,(x) and xG.(x) are fundamental 
solutions of Eq. (5), from which, by use of Eqs. (2) and (4a), the stress coefficients 
follow immediately. 

5. Determination of the deflection functions. The deflection w(x) can be expressed 
in the form 


w(x) = R[wi(x) + we(x) + wa(x)], (18a) 


where w;(x) and w(x) arise from the complementary functions respectively, and 
®s3(x) arises from the particular integral. The calculation of w#3(x) presents no diffh- 
culty, since only elementary functions with known integrals are involved. Direct in 


tegration of Eq. (6) gives 





2PR(1 + a) - — 3c x ] (19) 
1— <x 1—<x 


W3( x) ff eatayas TERA — 30) + log. ; 

The construction of the deflection functions w;(x) and we(x) is considerably more diffi- 
cult, since it is necessary to evaluate integrals of the type [xG;(x)dx (¢=1, 2), which 
involves additional infinite series. For purposes of computation, a convenient proce- 
dure, that also has the advantage of being easy to check, is to use a combination of 
analytical and numerical methods. A prerequisite for this calculation is a fairly ex- 
tensive and accurate tabulation of the G’s. 

A straightforward step by step numerical integration process is seen to fail near 
the poles of the G’s, due to the presence of ordinary singularities in ‘the integrands. 
The procedure for constructing the functions w;(x) and we(x) in tabular form con- 
sists of removing these singularities analytically and integrating the resulting func- 
tions numerically. 

6. Removal of singularities from the integrands of /xG,(x)dx and [xG2(x)dx. Let 
us consider a “substracting off” function H,(x) which has the property that 
G;(x) —H;(x) is bounded uniformly, i.e., without finite jumps, throughout the in- 
terval of existence of G;(x). It is necessary that H1;(x) be continuous except for poles 
which are of the same order as, and coincide with, those of G:(x). This specification 
is not sufficient however, since at every point of the interval the difference 
G(x) —H,(x) is finite, which requires the principal parts of Gi(x) and H,(x) to be 
identical. The principal parts of Gi(x) at zero and unity are readily obtainable from 
Eqs. (11) and (17) together with the relation G,(x)=G,(i-—x). Since G2(x) 
= —G,(1—x) and the F’s and G's are linearly dependent, the corresponding principal 


parts of G(x) may be found by the same process. 
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The integral parts of Gi(x) and G2(x) can. be approximated by polynomials of 
low ceyree, which makes it possible to reduce the differences Gi(x)—Hi(x) and 
G2(x) —H2(x) to uniformly small values throughout the interval by an intelligent 
choice of oe “subtracting off” functions. Incidentally this process provides a con- 
venient check on the accuracy of the tabulated values of the G’s. After some manipu- 


lation a pair of suitable “subtracting off” functions were found to be 


H,(x) = C2[— 7 + ab — g(0) + Ci/C2 — 4(ab — 4)(ab — 6) log, x(1 — x) 


+ 1/x? + 1/(1 — 2)? — (ab — 6)/x — (ab — 6)/(1 — 2)], (20) 
H(x) = D2[(7 — ab — g(0) + Di/Ds)(1 — 2x) — 3(ab — 4)(ab — 6) log, {x/(1 — x)} 
+ 1/x? — 1/(1 — x)? — (ab — 6)/x + (ab — 6)/(1 — 2)]. (21) 


These functions have the added property that 


lim {Gi(x) — Hi(x)} = lim {G.(x) — H2(x)} = 0. 


z—0, z—1 z—0, z—1 


Integrals of the type {xG(x)dx now can be evaluated directly from the identity 


— w(x) = f sGunae = f xe) — H;(x) lax + f xtts(a)as, (¢ = 1,2). (22) 


The second integral on the right-hand side is expressible in terms of elementary func- 
tions, while the first one behaves like a polynomial which can be computed easily 
with any numerical integration formula having a suitably small remainder depending 
on the magnitude of the differences of x[G;(x) —H,(x)]. Evaluation of the second 
integral of Eq. (22) with o=.3 gives, with the constant of integration chosen so that 


w:(4) = we(4) =0 
ff smioae = (060,042,744 .082,589,7:* + .052,500,0% — 2.103,471 


1 
oS ome + log. x — 1.047,500,0 log, (1 — x) — .052,500,0x? log, x(1 — ot, (23) 
— > 4 


ff etna = 040,784,504 .080, 308, 8 — .060,231,6x%7 + 4.147,500x + 1.349,471,6 


1 1—x 
_ aaa + log, x + 1.0475,000,0 log, (1— x)+ .052,500,0x? log, —*1. (24) 
—2x x 


7. Deflection and stress coefficients. It is convenient to state the actual deflec- 


tion in the form 


2R*(1 — 2 
w= noi Ahi J [. Aw, + Bwe + (P/ho)ws + Cc], (25) 


i Eho 
where, from Eqs. (18a) and (19), 


1 - — 3e 41 x =] 
Ws = ———— | ——— em , 
=I —oOLi-—s ee 


and w,; and we are non-dimensional functions of x defined by Eq. (22). wi, we, ws 
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are deflection coefficients, and are tabulated in Table 2. The constants A, B, and C 
are seen to have the dimensions of stress. 

The bending stresses can be stated in a form entirely analogous to that obtained 
for the rotating disc problem. With the aid of Eqs. (4a), (2), and (25), we have 


S, = Apit Bp2+ (P/h) ps, (26) 
Si = Aq + Baz + (P/K)4s, 
where 
pi = (1 — x) [xdG,/dx + (1+ 0)G,], 
gi = (1 — x)[oxdG,/dx + (1 + 0)Gi]; 


a similar pair of relations apply to p2 and g2; when ¢=.3, ps3 and g3 can be computed 
directly from the formulas 








1 — x [.63 + 2.27% — .7x? = 7 
ps = — 4.547 284 ——— pa | 
x (1 — x) x (27) 
1— «72.1 — 2.144 + .7x%? = .7 
gs = — 4.547 284 —|. 
x (1 — x)3 x 


Pi, Pe, Ps, G1,» Je, Ys are the stress coefficients, and are tabulated in Tables 3 and 4. 

The tables of coefficients are especially convenient for approximating a plate of 
variable thickness with a system of conical profiles. Calculations in this type of prob- 
lem show that it is necessary that the coefficients be accurate to six significant figures 
in order to obtain four significant figures in the final results. Consequently the tables 
have been calculated accurately to five parts in two million. Their general usefulness 
can be extended considerably with the aid of an auxiliary table of interpolation 
coefficients. It was found that such a table based on Bessel’s central difference 
formula® for six ordinates gives interpolated values of the coefficients as accurately 
as the tabulated ones, except near the ends of the table where the values are seldom 
used. In such cases a knowledge of the singularities of the tabulated functions indi- 
cates the necessary procedure for applying an interpolation formula. 


ILLUSTRATIVE EXAMPLE 


Stress distributions in a steel valve of constant weight and various proportions 
were estimated by an approximate method based on thin conical disc stress coeffi- 
cients tabulated for both the lateral bending and rotating cases. The valve head is 
represented by a system of truncated conical shells of variable thickness, whose apex 
angles are nearly 180° as shown in Fig. 1. The angle of the seat determines the direc- 
tion of the reaction which imposes two independent stress systems on the valve head. 
An approximation to these stresses can be made on the assumption that the mem- 
brane and bending stresses correspond to those in an equivalent system of conical 
discs. This assumption is admissible, since it has been demonstrated? for conical shells 
of constant thickness, that the stress distribution has the same character as that in a 





® J. B. Scarborough, Numerical mathematical analysis, The Johns Hopkins Press, Baltimore, Md., 
1930, p. 64. 
° Ref. 1, p. 477. 

















212 K. E. BISSHOPP [Vol. II, No. 3 


circular plate whenever the apex angle of the shell is between 168° and 180°. The 
loads on the composite disc shown in Fig. 2 are determined by resolving the valve 
seat reaction (of which the axial force is a component) into two perpendicular com- 
ponents, one of which produces pure compression on a section normal to the middle 




















Fic. 1. Half section of valve. 


surface at the periphery, and: the other of which produces pure bending.* The 
peripheral forces per unit length are proportional to the resultant vertical force 
acting on the valve, so that the force resolution in Figs. 1 and 2 has been made in 
terms of the axial force P, which is considered as a concentrated load, such as would 


PSin(@+x) 


_ Cos8 


—_— 
Cos 















Fic. 2. Half section of composite disc. 


be imposed on the valve by impact against its seat. It is safe to assume a concentrated 
axial load since the impact forces are proportional to the total valve weight, of which 
approximately 50% is in the stem. 

The next step in the calculation is to represent the valve head by a system of 
equivalent conical discs in the usual manner. The tabular solution for the bending 
stresses is obtained from the calculation procedure described in Ref. (2), except that 


* Variation in the slopes of corresponding generators of middle surfaces belonging to the conical discs 
of the equivalent system is not considered. 
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the p’s and q’s now refer to lateral bending coefficients and pRiw? is replaced by 
P/h.,. The solution for the membrane stresses is unchanged with the exception that 
the coefficients of p; and gs are zero, which corresponds to a static stress distribution 
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in a rotating disc. If S, and 5, refer to the corresponding membrane stresses respec- 
tively, then the appropriate boundary conditions are: at the boundary between valve 
stem and head, S,=oS,, 5,=S,; at the periphery S,=0, S, assigned. 


2.0- 


STRESS FACTOR S/P 
° a 
, _ 


w 
T 














1 
RADIUS IN INCHES 
Fic. 4. 


The dimensions of the valve head and the results of the stress calculations are 
shown in Figs. 3 and 4 respectively. 








TABLE 1.*—Fundamental solutions of hypergeometric equation, ¢ =.30. 
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x G(x) =Gi—x) | Gi(x)=—Gi(1—x) | Gi(x)=—Gi(l—x) | Gi (x) =Gz (1—x) 
.00 re) | — 0 0 —o 
01 613.258 —121,346.7 | 416.297 —82,426.1 
.02 156.631,5 | —15,325.94 | 106.126,0 —10,410.62 
03 71.137.6 —4'587.63 | 48.050,0 ~3'116.55 
.04 | 40.900,5 —1,955.005 27 .507,5 —1,328.327 
.05 | 26.762,0 —1,010.946 | 17.900 ,04 —687 .077 
.06 19.005 ,02 —590.775 | 12.627, 16 —401 .683 
07 14.282 ,02 —375.613 9.415 ,02 —255.543 
.08 11.187 ,24 —254.002 | 7.308 ,76 —172.948,9 
.09 9.045 ,62 —180.034,7 | 5.849, 83 —122.717,7 
.10 | 7.499,77 —132.421,0 4.795 ,46 —90.387,7 
11 6.345, 89 —100.354,8 4.007 ,22 —68.619,0 
12 5.460,74 —77.948,5 3.401 ,39 —53.412,3 
.13 | 4.766, 18 —61.804,5 2.924, 89 —42.459,9 
14 4.210,70 —49.867,4 2.542,72 —34.365,6 
15 3.759,18 —40.842,7 2.231,01 —28.250,3 
.16 3.386 ,98 —33.887,3 1.973 ,046 —23.541,2 
17 3.076,44 —28.435,8 1.756,797 —19.854,16 
.18 | 2.814,56 —24.098,7 1.573 ,446 —16.924,95 
.19 | 2.591 ,67 —20.602,3 1.416 ,397 —14.567,50 
.20 2.400 ,38 —17.749 , 86 1.280,639 —12.648,25 
21 2.235 ,01 —15.397,53 1.162,300 —11.069,57 
e 2.091,12 —13.438,55 1.058, 356 —9.758,98 
.23 1.965, 196 —11.792,44 .966,410 —8.661,87 
.24 1,854,433 —10.397,74 884,547 —7.736,55 
25 1.756,562 —9.206 ,96 811,215 | —6.950,83 
.26 1.669,738 —8.182,97 745,149 | —6.279,56 
27 1.592 ,445 —7.296,48 | .685 , 309 —5.702,90 
.28 1.523 ,429 —6.524,11 630 ,829 —5.205 ,08 
.29 1.461,645 —5.847,10 .580 ,987 —4.773,43 
| | ————— 
30 1.406,220 —5.250,21 535,174 —4.397,71 
31 1.356 ,416 | —4,721,02 492,874 —4.069,55 
32 1.311,608 —4.249,27 453,647 —3.782,13 
.33 1.271,267 —3.826,49 417,115 —3.529,80 
34 1.234,939 —3.445,58 382,950 —3.307 ,88 
.35 1.202, 236 —3.100,60 .350 , 869 | —3.112,48 
36 1.172,824 —2.786,52 | 320,623 —2.940,33 
.37 1.146,417 —2.499,05 .291 ,994 —2.788,70 
.38 1.122, 767 —2.234,53 264,788 —2.655,29 
.39 1.101,661 —1.989,781 .238 , 834 —2.536,15 
.40 1.082,915 —1.762,067 | .213 ,976 | —2.435 ,65 
41 1.066,371 —1.548 ,980 .190,076,6 | —2.346,40 
42 1.051, 893 —1.348 ,397 .167 ,007 ,9 —2.269,25 
43 1.039 ,367 —1.158,434 144 ,654,6 ~2.203 ,20 
44 1.028 ,695 — .977 ,394 .122,909,6 | —2.147 ,46 
45 1.019,795 — .803,742 | .101,673,3 | —2.101,34 
46 1.012,600 | — .636 ,068 .080,852,4 | —2.064,31 
47 1.007 ,058 — 473,063 | 060,358 ,2 —2.035,94 
.48 1.003 ,128 — 313,494 | 040,105 ,8 —2.015,90 
49 1.000, 780 — .156,186,3 | .020 013 ,2 —2.003 ,96 
.50 1.000 ,000 0 | 0 —2.000,00 





* The tables were compiled with the aid of the staff of the Calculation Deparisaent of Fairbanks 
Morse & Co., Beloit, Wis., to whom acknowledgement hereby is made. 


215 





















































1944] BENDING OF CONICAL DISCS 
TABLE 2.—Deflection coefficients for lateral bending of conical discs, ¢ =.30. 
r/R Ww We | W3 r/R W We Ws 
.00 x x | -* .50 0 0 10.00402 
01 .341490 .1727156 |—15.84278 51 — .00505132} .0000506834| 10.39010 
.02 . 298578 .1436068 |—12.59312 .52 — .01021073| ,000205606 10 .78484 
.03 .272917 .1262435 |—10.65009 .53 — .01548680; .000469405 11.18891 
.04 254305 . 1136970 —9.24108 .54 — .0208887 | .000847188 | 11.60306 
.05 . 239546 . 1037962 —8.12392 So — .0264263 | .001344576 12 .02809 
.06 .227215 .0955742 —7.19072 .56 — .0321101 | .001967757 12 .46484 
.07 .216551 .0885157 — 6.38390 oe — .0379516 | .00272355 12.91425 
.08 .207102 .0823131 —5.66908 .58 — .0439632 | .00361945 13 .37730 
.09 . 1985728 .0767682 —5.02404 .59 — .0501582 | .00466373 13 .85508 
.10 . 1907641 .0717460 | —4.43361 .60 | —.0565514 | .00586555 14.34880 
as . 1835327 .0671501 —3.88695 61 — .0631588 | .00723498 14.85973 
a E .1767727 .0629098 —3.37604 .62 — .0699978 | .00878323 15 .38931 
13 . 1704034 .0589712 —2.89476 .63 — .0770879 | .01052271 15 .93909 
14 | .1643621 .0552928 —2.43835 .64 — .0844501 | .01246720 16.51082 
15 | .1585986 .0518418 — 2.00300 65 — .0921083 | .01463208 17.10641 
16 | .1530724 .0485919 —1.585657 .66 — .1000884 | .01703454 17 .72799 
17 .1477501 .0455217 —1.183777 .67 — .1084197 | .01969381 18.37793 
.18 . 1426040 .0426135 — .795248 .68 — .1171349 | .0226315 19 .05890 
19 | .1376105 .0398528 — .418284 .69 — .1262706 | .0258720 19 .77389 
20 | .1327494 .0372275 — .0513587| .70 — .1358681 | .0294430 20.5263 
21 1280035 .0347271 .306849 71 — .1459743 | .0333756 21.3199 
‘ae .1233577 .0323431 .657492 me: — .1566424 | .0377057 22.1591 
‘an . 1187984 .0300681 1.001586 3 — .1679333 | .0424741 23 .0488 
24 . 1143140 .0278957 1.340034 74 — .1799171 | .0477280 23.9948 
.25 | .1098938 .0258206 1.673648 .75 — .1926744 | .0535219 25 .0038 
.26 | .1055282 .0238382 2.00316 .76 — .206299 .0599190 26.0833 
.27 . 1012084 .0219446 2.32924 ae — .220901 .0669937 27.2424 
.28 | .0969264 .0201366 2.65250 .78 — .236610 .0748334 28.4918 
.29 | .0926748 .01841116 2.97352 .79 — .253578 .0835418 29.8439 
.30 .0884466 .01676613 3.29283 .80 — .271990 .0932430 31.3139 
31 .0842354 .01519954 3.61092 1 — .292065 . 1040868 32.9200 
.32 .0800348 .01370985 3.92829 .82 — .314074 . 1162557 34.6842 
33 .0758391 .01229585 4.24537 .83 — .338347 .1299744 36.6339 
34 .0716426 .01095665 4.56262 . 84 — .365296 | .1455227 38 .8030 
35 .0674398 .00969164 4.88046 .85 — .395443 . 1632539 41.2345 
.36 .0632252 .00850050 5.19930 .86 — .429454 . 1836207 43 .9833 
an .0589937 .00738318 5.51956 .87 — .468197 .207213 , 47.1212 
.38 .0547400 .00633986 5.84165 .88 — .512828 . 234815 50.7436 
.39 .0504589 .00537098 6.16596 .89 | —.564915 267491 54.9800 
.40 .0461453 .00447722 6.49292 .90 — .626646 | .306728 60.0115 
41 .0417938 .00365951 | 6.82293 91 —.701169 | .354663 66.0987 
.42 | .0373993 .00291902 7.15642 .92 — .793190 .414495 73.6312 
.43 | .0329561 .00225716 7.49380 .93 — .910079 .491227 83.2197 
44 .0284588 .001675592| 7.83553 .94 |—1.064072 .593171 95 .8789 
45 | .0239016 .001176243| 8.18206 95 |—1.277112 . 735232 113.4294 
.46 | .01927865 .000761311) 8.53386 -96 |—1.592926 .947112 139.5018 
.47 | .01458362 .000433277 8.89143 97 |—2.11315 1.297850 182.5406 
.48 | .00981010 .000194922/6 9.25528 .98 |—3.14154 1.993762 267 .798 
49 .00495128 .000049349/6 9.62595 .99 |—6.19036 /|4.06209 521.097 
0 10.00402 | 1.00 —« ~ “ 
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TABLE 3,——Stress coefficients for lateral bending of conical discs, o =.30. 

' o/R | bi -p | pb |r/R|  p -p | bs 
.00 —« | 20 r | .50 .650000 | —_ .500000 —51.4753 
.01 | —412.069 280.244 '31209.9 | .51 676528 513539 —53.6293 
.02 |—100.8399 68.8436 7638.76 | .52 . 704200 528194 —55.8788 
.03 | —43.7956 30. 1006 3318.31 .53 733152 544032 —58.2347 
04 | —24.0284 16.67837 1821.085 .54 763534 561125 —60.7094 
.05 | —14.96886 | 10.52962 1134.831 .55 795506 .579561 — 63.3159 
.06 | —10.09557 | 7.22454 | 765.640 .56 829243 .599438 — 66.0685 
07 | —7.18545 5.25310 545.144 57 864938 .620867 — 68.9831 
.08 | —5.31459 | 3.98776 403.367 .58 .902803 643975 —72.0770 
.09 | —4.04388 3.13023 307.050 | .59 .943074 668906 —75.3697 
10} —3.14316 | 2.52420 238.760 | .60 .986012 695824 — 77.8828 
"11 | —2.48255 | 2.08145 188.6601 | .61 1.031911 .724915 — 82.6402 
12 | -—1.984277 | 1.749150 150.8584 | .62 1.081101 756391 — 86.6694 
13 | —1.599538 | 1.494168 121.6581 | .63] 1.133955 . 790495 ~91.0010 
14] —1.296465 | 1.294864 98.6451 | .64 1.190896 .827503 —95 .6698 
15 | —1.053553 | 1.136640 80.1902 | :65 | 1.252404 .867734 | —100.7157 
16 —.855874 | 1.009365 65.1626 | .66| 1.319031 911553 — 106.1840 
17} —.692812 |  .905838 52.7577 | .67| 1.391409 959381 112.1271 
18 | —.556649 .820830 42.3909 | .68 1.470270 1.011709 — 118.6053 
19 | — .441674 . 750472 33.6293 | .69 1.556461 1.069105 — 125.6887 
.20| —.343586 | .691855 26.1469 | .70 1.650971 1.132236 — 133.4592 
21] —.259098 | .642759 19.69457 | .71 1.754958 1.201882 — 142.0124 
.22 —.1856618 |  .601468 14.07923 | .72 1.869789 1.278967 — 151.4616 
.23 | —.1212805 |  .566640 9.14940 | .73 1.997083 1.364585 — 161.9406 
.24| —.0643692 | .537215 4.78487 | .74 2.13878 1.460047 — 173.6096 

| | | 
.25 | — .01365723| 512346 39247| .75 2.29719 1.566925 — 186.6610 
.26 .0318840 491353 —2.61552 | .76 2.47513 1.687126 —201.327 
.27 0730951 | 473684 —5.79325 | .77 2.67604 | 1.822974 — 217.893 
.28 . 1106690 .458889 —8.69658 | .78 2.90412 1.977330 — 236.707 
.29 1451814 .446599 | —11.36928 | .79 3.16461 | 2.15375 — 258.203 
.30 1771155 | .436510 | —13.84808 | .80 3.46408 | 2.35669 — 282.926 
31 .206879 | .428369 | —16.16404 | .81 3.81084 2.59179 —311.566 
32 .234821 | .421968 | —18.34365 | .82 4.21558 2.86631 — 345.007 
.33 .261238 | .417132 | —20.4097 | .83 4.69218 3.18967 — 384.402 
.34 286390 | .413718 | —22.3820 | .84 5.25895 3.57432 —431.270 
| 

35 | 310503 | .411605 | —24.2778 | .85 5.94048 4.03696 —487.652 
.36 | .333776 | .410694 | —26.1125 | .86 6.77038 4.60040 —556.337 
.37 | .356387 | .410903 | —27.8996 | .87 7.79558 5.29653 —641.223 
.38 | 378496 | .412167 | —29.6516 | .88 9.08323 6.17096 —747.885 
.39 | 400248 | .414431 | —31.3796 | .89| 10.73220 7.29083 — 884.537 
40 | 421777 | .417655 | —33.0941 | .90| 12.89286 8.75831 — 1063 .672 
41 | .443208 | .421806 | —34.8048 | .91| 15.80318 10.73501 — 1305 .066 
42 | ‘464658 | .426865 | —36.5208 | .92| 19.85801 13.48915 —1641.541 
43 486239 | .432816 | —38.2510 | .93 | 25.7521 17.49263 — 2130.86 
44 | .508060 .439656 | —40.0040 | .94| 34.8021 23.6398 — 2882.53 
45 | 530227 | .447386 | —41.7881 | .95 | 49.7595 33.7997 —4125.41 
.46 | 552846 | .456017 | —43.6117 | .96| 77.1990 52.4381 — 6406.56 
47 | .576023 | .465566 | —45.4835 | .97 | 136.2744 92.5655 —11320.03 
.48 599866 |  .476057 | —47.4119 | .98 | 304.461 206 . 807 —25315.3 
49 .624487 | .487522 | —49.4061 | .99 |1209.305 821.430 — 100647 .9 
-50 | .650000 | .500000 | —51.4753 [1.00 0 20 | — © 
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TABLE 4.—Stress coefficients for lateral bending of conical discs, ¢ =.30, 



























































r/R | q1 qQ: —s r/R qQ | q2 —Qs 

00 | 00 « “ 50 .650000 | —.1500000/ 50.2020 
01 |428.863 290.968 32477.1 51 .649206 | —.1629858} 50.2128 
02 |109.4320 73.9901 8285 .68 .52 .649426 | —.1759765| 50.3035 
.03 | 49.6545 33.3836 3758.74 53 .650664 | —.1890247| 50.4748 
.04 | 28.5222 19.02706 2158.48 54 652935 — .202182 50.7282 

| 

05 | 18.64509 12.31570 1410.590 55 656258 — .215504 51.0654 
06 | 13.22822 8.63391 1000 .466 56 660662 — .229044 51.4891 
.07 | 9.93124 6.39200 750.873 .57 666186 — .242863 52.0026 
08 | 7.77158 4.92257 587.402 .58 .672875 — .257023 52.6097 
09 | 6.27752 3.90518 474.331 .59 .680785 — .271589 53.3154 
.10 5.19936 | 3.17022 392.751 .60 689984 — .286635 54.1253 
ll 4.39477 | 2.62101 331.883 61 700553 — .302237 55.0462 
12 | 3.77767 | 2.19909 285.210 62 712583 — .318481 56.0859 
13 3.29352 | 1.867385 248.602 .63 726185 — .335463 57.2535 
14 2.90636 | 1.601472 219.336 .64 741486 — .353287 58.5598 
15 | 2.59166 | 1.384697 195.5559 .65 . 758634 — .372072 60.0169 
16 | 2.33225 | 1.205387 | 175.9608 | .66 .777800 — .391950 | 61.6393 
17° | 2.11579 | 1.055157 159.6167 .67 799185 — .413074 63.4435 
18 | 1.933233 .927856 145 .8386 .68 823021 — .435615 65.4490 
19 | 1.777816 | .818885 134.1146 .69 849583 — .459771 67.6783 
.20 1.644399 .724748 124.0556 .70 879189 — .485773 70.1581 
21 1.529018 .642750 115.3617 71 912216 — .513887 72.9195 
22 1.428577 570780 107.7985 .72 949106 — .544425 75.9993 
23 1.340629 | .507172 101.1808 .73 .990389 — .577756 79.4413 
24 1.263215 |  .450588 95 .3605 74 1.036693 — .614316 83.2977 
25 1.194757 .399950 90.2181 .75 | 1.088774 — .654629 87.6312 
.26 1.133967 354377 85.6563 .76 1.147548 — .699322 92.5176 
27 1.079790 313145 81.5951 77 1.214126 — .749161 98.0492 
.28 1.031351 .275652 77.9685 .78 1.289876 — .805081 | 104.3394 
.29 | .987924 . 241396 74.7215 .79 1.376492 — .868241 | 111.5287 
30 |  .948897 . 209953 71.8079 .80 1.476092 — .940082 | 119.7926 
31 913757 .1809649 69.1890 81 1.591352 —1.022432 | 129.3533 
32 .882070 1541266 66.8319 .82 1.725701 —1.117624 | 140.4952 
33 853463 .1291751 64.7086 83 1.883578 —1.228709 | 153.5868 
34 827622 1058843 62.7951 84 2.07083 —1.359573 | 169.1134 
35 804274 0840574 61.0712 .85 2.29527 —1.515621 | 187.7236 
36 | .783186 0635226 59.5192 .86 2.56756 —1.704062 | 210.302 
37. | .764157 0441290 58.1239 .87 2.90251 —1.934972 | 238.080 
38 | .747014 0257434 56.8723 .88 3.32128 —2.22272 272.814 
.39 | .731606 00824728} 55.7532 .89 3.85488 —2.58838 317.080 
.40 .717805 | —.00846535| 54.7569 .90 4.55034 —3.06387 374.786 
.41 | .705497 | —.0244898 53.8749 91 5.48179 —3.69960 452.095 
42 | .694587 | —.0399126 53.1001 .92 6.77183 —4.57882 559.197 
43 | .684992 | —.0548125 52.4264 .93 8.63538 —5.84753 713.963 
.44 | .676641 | —.0692620 51.8486 .94 | 11.47831 —7.78139 950.147 
45 | .669475 | —.0833284 51.3624 .95 | 16.14551 —10.95435 | 1338.035 
.46 | .663446 | —.0970743 50.9640 .96 | 24.6485 —16.73272 | 2045.00 
47 |  .658511 | —.1105590; 50.6507 .97 | 42.8244 —29.0814 3556.83 

.48 | .654640 | —.1238390 50.4203 .98 | 94.1889 —63.9738 7831.14 

49 | .651808 | —.1369684 50.2710 .99 |368.372 —250.217 30658.4 

.50 | .650000 | —.1500000 50.2020 | 1.00 20 —« x 

















KRON’S METHOD OF SUBSPACES* 


BY 
BANESH HOFFMANN 
Queens College, Flushing, N. Y.** 


Introduction. Gabriel Kron has introduced new and powerful methods of applying 
tensor analysis to complicated engineering problems, presenting his major contribu- 
tions in the field of electrical engineering. The manner of presentation, and the rarity 
of a simultaneous knowledge of the hitherto almost unrelated subjects of electrical 
engineering and tensor analysis, have unfortunately served to limit his audience. 
Recent experimental confirmation of some of his investigations dealing with equiva- 
lent circuits, however, has attracted the serious attention of a wider engineering 
following. 

In view of the growing importance of the whole subject, and of the controversy 
which has surrounded it, it has seemed desirable to present some particular aspect 
of Kron’s work in a form which may appeal to a less highly specialized audience. To 
avoid complications as far as possible, the present paper must ignore such important 
topics as electrical networks, electrical machines, and equivalent circuits. It confines 
itself to purely dynamical problems, and to that particular idea of Kron’s which may 
be called the method of subspaces. 

Much discussion has arisen over Kron’s claim that he uses tensor analysis. It is 
the considered opinion of the present writer that Kron does indeed make a full and 
proper use of tensor analysis. Possibly the belief that Kron employs only matrices 
may have arisen from the fact that, in order to present his actual mathematical 
procedure in a form that may be understood and used by those not familiar with the 
intricacies of the tensor calculus, he often presents this procedure in matrix form. 
However, he is always careful to point out that the underlying concepts are wholly 
tensorial in character. In the present paper the method of subspaces will first be pre- 
sented in terms of a simple example, and in purely matrix form, merely as a set of rules 
of procedure, the essentially tensorial significance of the procedure being discussed 
only after the actual procedure has been brought before the reader. The theoretical 
discussion will then be followed by three simple, related examples illustrative of 
various aspects of the method. 

Scope of the method of subspaces. Let us consider a system, which may be dy- 
namical, electrodynamical, or otherwise, containing several standard parts, such as a 
fly-wheel, a governor, a pair of synchronous machines, a system of levers, etc. The 
equations of performance of the individual parts are usually well known, but the 
equations of performance of the complex whole will depend on the manner in which 
they are interconnected. Usually it is extremely difficult to trace out the full influ- 
ence of each interconnection in setting up the equations of performance. The method 


* Received Feb. 9, 1944. 
** On leave of absence; now at Federal Telephone and Radio Laboratories, New York City. 
+ Purely dynamical examples of the method of subspaces have been given by Kron in an unpub- 


ished manuscript. 
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of subspaces suggested by Kron yields the equations of performance by a routine and 
quite straightforward manipulation of the known equations of performance of the 
constituent parts of the system, a brief inspection sufficing to yield all needed informa- 
tion as to the manner of interconnection. 

Simple dynamical example. To illustrate the actual mathematical procedure in 
its simplest form, we shall first consider a quite trivial dynamical problem. Naturally 
it will not reveal the power and economy of the method any more than it would the 
peculiar virtues of, say, the Hamilton-Jacobi equation, were that applied to it. But 
it will serve to bring the routine mathematical procedure before us without unneces- 
sary distraction from complexities which are merely incidental. 

Let us consider the dynamical system S consisting of three particles free to move 
on a line, the masses being mm, !%2, m3, the coordinates! x’, x*, x*, and the forces fi, fo, fs. 
The equations of motion are 


fi = mx, fo = mx, Ss = msx*. (1) 


Let us consider now the new dynamical system 5 which arises when the particles 
2 and 3 are made to coalesce. Its equations of motion may be written down at once: 


fi = mz, Sot fs = (m2 + ms)z*. (2) 


Let us suppose, though, that it had been a highly complex system of interconnected 
simple parts. We would then welcome a routine method of obtaining (2) from (1) 
which required no detailed thought and avoided constant preoccupation with the 
effects of the interconnections. The method of subspaces would be applied to the 
present problem in the following routine manner: 

The first step is to write equations (1) of system S in matrix form, 


F= MX, (3) 
1.€., 
my, 0 0 x! ti 
0 me 0 #? = te i (4) 
0 0 m3 x* fs 


where, it will be noted, the masses form a square array rather than a single row or 


column such as one might at first expect. 
Next, a relationship is set up between the coordinates x!, x*, x* of S and the co- 


ordinates #!, £2 of 5. This relationship can be taken to be 


xz) = gi, x? = £%, x* = #? (not #°). (5) 


From this is obtained the matrix C defined by 


X = CX. . (6) 
It is 
1 0 
Cal @ i iL (7) 
0 1 


1 We use the tensor practice of placing the indices of contravariant quantities above the symbol. 
x!, x2, x8 do not stand for x, x-squared, x-cubed. 
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If we denote the transposed matrix by C“, then the new forces hi, fe are given by the 
routine matrix multiplication 


1 
- 7p 100 J fi 
F = CF = fe |= / ; (8) 
01 1 fotfs 
fs 
The new masses are given by routine matrix multiplications as follows: 
m O 0 1 0 
_ 1 0 0 
M =C"MC = 0 m, O 0 1 
01 1 
0 O ms 0 1 
1 0 
m, O 0 mM, 0 
= 01i{= b (9) 
0 m2. ms 01 O m+ ms 
) 


Finally the new equations of motion are 
F = MX, (10) 


PF emlleriace! 
0 Me + ms 52 7 he +t fs 


This yields the two equations 


2, 


mk = fy, (m2 + ms)#* = fo + fa, 


which are equivalent to (2) above. 

The tensor form of the problem. Using Latin indices for the range 1, 2, 3, and 
Greek for the range 1, 2,.we may write the various expressions and equations above 
in the familiar index notation of the tensor calculus. 

The equations of motion of S may be written? 





fa = Marx’, (1’) 
the matrix C may be written 
, ee (79 
c= , 
O#* 


and the relations between the forces, etc., in S and S may be put in the form 


; a (8’) - Ox* dx? 9’) Ox? 6’) 
eo” ay , Map = —— Mab, ’ xe = £%. pe 
oe aKa o ane ane” age 











Also the equations of motion of S are 
fa = Maph. (10’) 
Provided we interpret the time derivatives as absolute derivatives, or alterna- 


2 According to the summation convention of the index notation, a repeated index in a single ex- 
pression indicates summation over its whole range of values. 
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tively, as is permissible in the present simple case, provided we avoid non-linear 
coordinate transformations, the above have the form of tensor equations and ten- 
sor transformations. The mere fact that the work can be expressed in the tensor 
notation does not, in itself, imply tensorial character. The essential criterion of tensor 
character is the tensor law of transformation. In view of (7’), equations (8’), (9’), 
and (6’) are tensor transformations, and equation (10’) is the result of transforming 
the tensor equation (1’). The fact that C is singular does not destroy the tensor 
character of the transformations; its significance will appear shortly. 

The tensor theory. Since tensor equations have objective significance we may 
look for a geometrical picture of the process described above. Naturally this will be 
sought in configuration space. The basic tensorial and geometrical significance of 
Lagrangean dynamics being quite familiar, the corresponding significance of Kron’s 
method may be explained here quite briefly. 

As is well known, in Lagrangean dynamics the motion of the dynamical system S 
is represented by the motion of a point in a three dimensional configuration space, K, 
having x* as coordinates and m,,» as metrical tensor. When particles 2 and 3 coalesce, 
the system S loses one degree of freedom and becomes the system S. Thus the 
trajectory of 5 belongs to a two dimensional configuration space, K, which is in fact 
a subspace of K, for it is defined by a relation (in more general cases, by a set of 
relations) between the coordinates of K. Specifically, the subspace here is defined by 


the relation 


In parametric form this subspace is given by the relations (5) above, which represent 
the three coordinates x* as functions of the two variables #*. (Compare with the rela- 
tions x =cos @ cos ¢, y=cos 6 sin yg, z=sin @ which express the Cartesian coordinates 
x, y, z as functions of the two parameters 0, g. This defines the two dimensional sub- 
space of ordinary three dimensional space constituting the surface of a unit sphere.) 

By the well known theory of subspaces, the projections of the covariant tensors 
fa, Ma, are given by the singular transformations (8’), (9’). Thus the equations of 
motion of 5 are the projection on K of the equations of motion of S. And, since the 
initial conditions of S and S coincide in K, the trajectory of 5 is the projection on K 
of the trajectory of S. 

There are two ways of viewing the relationship between the systems S and S: 

(a). We may regard S as the same physical system as S, the forces between 
particles 2 and 3 which keep them together being included explicitly in the force 
vector f,. The trajectory of S in K is then identical with that of S in K. 

(8). We may regard S as a different physical system from S inasmuch as 
particles 2 and 3 are not united in S. The forces in S are the same as those in S 
except for those forces in S which tend to separate particles 2 and 3. These 
latter forces have components in K which are normal to the subspace K, and thus 
have zero projection on K. 

Both viewpoints are of significance, and more will be said about them later. 

A formal proof. A formal proof will now be given that the transition to a sub- 
space and the tensor transformations that go with it, are justified in the general case. 
This proof will thus also justify the general procedure given by Kron, of which the 
above example was a particular illustration. 
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The proof will be made brief by basing it on certain standard results in dynamics 
and tensor analysis. 

Let us consider a rigorous proof of the validity of the Lagrangean equations, such 
as is given, for instance, in Whittaker’s Analytical Dynamics, third edition, starting on 
page 34. We are concerned here with the case in which ¢ does not enter explicitly, 
and we shall use the tensor notation. The dynamical system under consideration in 
Whittaker has m degrees of freedom, and n generalized coordinates. Let us denote 
the latter by #%, using the Greek indices a, 8, y for the range 1 to m. We denote the 
number of individual particles in the system by N/3, so that their combined coordi- 
nates number WN and we denote these N coordinates by x*, using \, wu, v for the range 
1 to N. In this notation, the proof of the Lagrangian equations has the following out- 
line: 

The equations of motion of the NV/3 individual particles, in a self-explanatory 
notation, have the form ‘s ae 
My, X" = p. (11) 
These NV equations are not independent, since the N coordinates x* are related, as 
are the N forces f,. The relations between the coordinates X* are defined by equations 


DP = X( #4 (12) 


which express them in terms of the ” generalized coordinates of the system. (These 
correspond to Whittaker’s equations x;=f;(q:, g2,---, Qn, t), etc., with ¢ omitted.) 
From (12) may be computed the quantities 0x*/0z%*. The equations of motion (11) 
of the individual particles are multiplied individually by such quantities and then 
added in groups, the process being precisely that described by the equation 


Se . fF. 
My, x" = — fy, (13) 





where the summation convention is employed, as usual. After some manipulation, 
the left-hand side is then reduced to the standard Lagrangean form, and the right- 
hand side is interpreted as a set of generalized forces in the familiar manner. 

Later, on page 39, Whittaker gives an explicit form of the Lagrangean equations 
for the case in which ¢ does not enter explicitly. This reveals that the left-hand side 
has the form of a covariant derivative with respect to #iag as metrical tensor. The 
equations, in fact, may be written 


Miah? ,427 = fa. (14) 


Since, in the original form (11), the coordinates were Cartesian for each individual 
particle, the ordinary derivatives there coincided with the covariant derivatives; thus 


(11) may be written as le os 
My," Xe” = Ty (15) 


the subscript preceded by a comma here denoting the covariant derivative with re- 


spect to m), as metrical tensor. 
It will be observed that the initial step, represented by the equation (13), is a 
transition to a subspace, complete with a singular transformation matrix 0x*/dz* 
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of the type which, for some reason, excites the critics of Kron when it is used by him. 
The relations between m7, and #i,g, and between f, and f,, follow from the analysis 
and are the usual tensor transformations, with singular transformation matrix. 

Now Lagrange goes from an initial configuration space of N dimensions to a final 
subspace of » dimensions in one step. Kron’s idea is, in theory, simply to make this 
transition in more than one step, using subsidiary subspaces as resting places when 
the mathematics tends to become too complicated. For example, he would, in theory, 
go from the first space having coordinates x** to an intermediary subspace having 
coordinates x* and then to the final subspace having coordinates #* which is also a 
subspace of the intermediary subspace. In practice, of course, as is also the case in the 
Lagrangean method, the initial space having coordinates x is entirely neglected, 
having served its sole purpose in providing the theoretical basis for the equations and 


procedures actually used. 


Since’ 
Ox dx* Ox 


ax® af" axe 











b 


the transition in several steps will yield the same result as the transition in one step, 
the various tensors involved being transformed according to the standard tensor law. 

Thus the proof of Kron’s theory and procedure is a direct corollary of the proof 
of Lagrange’s equations and the tensor theory of subspaces. 

The proof given in Whittaker is based on viewpoint (a), since the forces f, be- 
tween the individual particles are regarded as the forces actually existing between 
them in the ultimate system. The trajectory in the N dimensional space is actually 
confined to the m dimensional subspace. 

It is important to note, though, that the proof is equally valid for viewpoint (8). 
For those forces which do no work do not contribute to the values of the generalized 
forces of the ultimate system. Since they do not affect the ultimate system, it is clear 
that, for the purpose of setting up the equations of that system, they may be omitted 
from the N basic equations of the individual particles. When these forces are ignored, 
however, the forces between the individual particles are very much changed, espe- 
cially in the case of inelastic bodies. The system of individual particles is then no 
longer physically equivalent to the ultimate system. It has a quite different motion 
(for instance, the particles of a rigid body here move in divergent directions) and its 
trajectory, in general, spans the whole N dimensional space. Nevertheless, according 
to the above reasoning, the projection of its trajectory on the m dimensional sub- 
space coincides with the trajectory of the ultimate system. 

Non-linear transformations. Since Kron has made the widest application of his 
method of subspaces to electrical networks and other electrodynamical problems in 
which the interconnection transformation is very often linear, the impression has 
sometimes arisen that the method is applicable only to situations in which this 
linearity is present. The following examples of the method involve non-linear trans- 
formations. They are simple enough so that the more usual methods of solution are 


3 When the transformations are non-singular, this is the basis of the important “group property” of 
tensor transformations. In the absence of inverses the term “group” is inappropriate here, but provided 
the succession is always to a subspace of the preceding space or subspace, as it always is here, the 
usual combination properties of tensor transformations are preserved, 
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hardly more complicated than those by Kron’s method, but this is inevitable in simple 
problems. The power of Kron’s method begins to be felt when the systems in question 
involve a larger number of interconnections, and of more complicated mechanisms 


than the simple rods of the examples below. 
Illustrative example I. We begin with 





Z 
1 a system consisting of two rods hinged 
together without friction, one rod being 
0 ™ r 
JY suspended by its free end from a fixed 
point O (Fig. 1). We denote the masses 
8 of the rods by m, me, their lengths by 


2a;1, 2a2, and their moments of inertia 
| about their centers of gravity, which we 


1A shall assume to coincide with their mid- 

mM, g iP F points, by J;, J». In addition to gravity, 
! : B let us consider a force F, not necessarily 

! conservative, which acts horizontally and 

mg to the right at the mid-point of the lower 


rod. The system has two degrees of free- 
dom, and we may take as the generalized 
coordinates the angles #0, g which the 
Fic. 1. System I. rods make with the vertical. 
The problem is to set up the equations 





of motion. From previous experience with Lagrangean dynamics, we may regard a 
single rod as a known system, in the sense that we can instantly write down its equa- 
tions of motion, or have already tabulated them for quick reference. The present 
system consists of two of these 
known systems interconnected. We 
therefore begin by considering the 
system consisting of the two rods O ~y 
not interconnected, the forces being 
the same as those acting externally 6 
on the original system.‘ Kron calls 
this system the primitive system. 
It is shown in Fig. 2, and has here 
four degrees of freedom. We may A! 
take the four generalized coordi- mg lo F 
nates to be the angles 0, ¢ above 
together with the coordinates y, z 

of the center of gravity of the lower , M28 
rod. We let Latin indices refer to 
the primitive system, and Greek 
to the actual system under discus- = 

; rm ea ge Fic. 2. Primitive system of system I. 
sion. For the primitive system the 

metrical tensor and the force vector can be written down at once. They are 


IN 


me | 








* In general one must include all forces that do work, including dissipative forces. 
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* fe fel 
@ | 7, +m,a2 0 | 0 | 0 
me = go | 0 I | 0 | 0 (16) 
y | 0 0 | ms | 0 
ra 0 0 oi 
and 
6 | — mga; sin @ | 
| 
¢ | 0 
fo = (17) 
y | F | 
2 | — M2g 





The restraint arising from the interconnection of the two rods imposes the following 
two conditions on the four coordinates x* of the primitive system: 


y = 2a, sin 6 + az sin ¢, z =— 2a, cos 6 — a: cos¢. (18) 


These two equations define the subspace of the configuration space of the primitive 
system to which the given system is confined. We may express them in the form of a 
transformation, that is to say, in parametric form, by writing 























6=6, e=%, y=2a,sin6+a:sing, z= — 2a,cos§—a,cos¢, (19) 
which is of the form 
x* = x*(2*), (20) 
The transformation matrix C2, or 0x*/0#+, is (in the form C”) 
x 
pe ° ¢ y a 
Ox? — = = 
= = 6 1 | 0 2a, cos 6 | 2a; sin @ |. (21) 
a 
rr) 0 | 1 a2 COS @ a2 sin @ 





Thus the metrical tensor for the given system, namely the projection Mas of ma, is 
given by 





a Ox* dx 
| Map = ee on 
i.e., 
Ih+mai 0 0 0 1 0 
_ 1 0 2a, cos 8 2a; sin 0 I, O 0 0 1 
ae i 1 ade:cos@ azsin -| 0 0 mm: O 2a,cos6 a:cosg 


0 0 O mz: 2a,sin6 assing 








226 BANESH HOFFMANN [Vol. II, No. 3 


1 0 
[” +m,a; 0 2aymz.cos 6 2a,m:2 sin 6 0 1 
- 0 Io  @oM2COSG aymz sin i cos6 a2:cos¢@ 
2a,sin6 a,sing 
7 [* + ma; + 4m2a?_ = 2mea 102 cos (8 — 9) (22) 
= 2m2a;a2 cos (8 — $) To + mea | - 


Likewise, the force vector is given by 


‘ Ox? 
oe = aze a) 
or ” 
— mga, sin 0 
>; _[1 0 2a, cos@ 2a; sin 8 0 
lo F 1 a.cos¢g asin “| F 
— Mg 
_ - (m; + 2mz)ga; sin 6 + 2a,;F cos 8 |. (23) 
a2F cos § — aemog sin G 


The kinetic energy function of the given system is 


T= 1 tig X*X8 

= 3(1, + ma? + 4mza?)b? + 2mara, cos (8 — BOG + 43(I2 + ma?) o. (24) 
The Lagrangean equations for the given system may now be written in the usual 
form.® They are, on dropping the bars over @ and ¢, but without simplification, 


d 
a {(1, + ma; + 4m,0?)6 + 2meaa2.cos (6 — 9) ¢ } — {— 2meaya, sin (0 — y)bo} 
= — (m, + 2m)ga, sin 0 + 2a,F cos 8, 


d 
= { 2ma,a2 cos (06 — ¢)6+ (I, + m2) ¢} — {2maaz sin (0— ¢)6¢} 
= a2F cos ¢ — d2meg sin ¢. 


Since no forces were introduced at the points A, A’ of the primitive system, the 
latter was physically different from the given system, for in the given system opposite 
forces acted at the hinge A. Thus we have been using viewpoint (8)..To make the 
two systems physically equivalent, it would be necessary to impose appropriate initial 
conditions on the primitive system and to introduce the proper opposite forces at 
A and A’ corresponding to the reactions at the hinge in the given system. This, 


5 Kron has suggested another method having advantages when the system and its interconnections 


are complicated. 
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however, would entail a knowledge of the reactions at A, and the advantage of the 
(8) viewpoint, which makes it the more appropriate viewpoint for the Kron method 
here, is that it enables one to proceed without bringing in the reactions at all, except 
indirectly insofar as they imply the equations of constraint. It is possible to use view- 
point (a) by introducing unknown opposite reactions at A and A’ in the primitive, 
denoting them by some symbol, say R and —R; they will automatically cancel when 
the transition is made to the subspace. 











To give some indication of the flexi- Z 
bility of the Kron method and its ability 2c -| 
to extract cumulative dividends from (8) ! y 
such calculations as may previously 
have been performed, we conclude with e 


a brief and sketchy discussion of two 
further systems. 

Illustrative example II. Let us con- 
sider the system illustrated in Fig. 3. mg 
It is the same as system I above except ' 
that the end of the lower rod is con- 
strained to move without friction on a 
fixed vertical line distant 2c from 0. 

System I has already been investi- 
gated. It is now a known system. In- ; 
stead, therefore, of taking the primitive Fic. 3. System II. 
of system II--to be the same as the 
primitive of system I, we may take it to be system I itself. 

The new constraint reduces the number of degrees of freedom to one and may be 
represented mathematically by the condition 


6|/> 





m8 





a; sinO+ asing =. (25) 


The subspace now can be written parametrically as 


6 = 6, y= sin f— E : (26) 


a2 


The transformation matrix C is given by 


1 ri) — a, cos 8 
c=| |. sai ! ; (27) 
N a V a? — (c — a; sin 6)? 








By implicit differentiation of (25) we may also obtain the useful relation 
a; cos 86 + asd cos g = 0. (28) 


The new metrical tensor #,, (which here has only one component) may be obtained 
by the usual transformation formula, or the new expression for T may be obtained 
directly from (24) by substituting for g and ¢ in terms of 9 and 6 by means of (26), 
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the latter method being the simpler in this particular problem. The generalized force is 


i a (m, + 2mz)ga, sin 6 + 2aF cos | 


a2F cos g — d2meg sin ¢ 


f. 


= — (m, + 2mz2)ga, sin 6 + 2a,F cos 8 + Aa2F cos o — Aagmeg sin ¢ 


aymee(c — a, sin 8) cos8 
1 2g ( 1 p) (29) 





— (m, + 2m2)ga; sin 6 + aiF cos 6 + 


oe eee , 
Vai — (¢ — a; sin 6)? 


the terms in F partially cancelling in view of (28). The equation of motion may now 
be written down in the usual manner. Solving it is another matter! 

Illustrative example III. The preceding example made use of system I in the role 
of a known system, and essentially dealt with the imposition of a constraint on that 
system. One may, however, join several known systems together by the Kron method, 
and this is, in fact, the procedure of principal importance in practical problems. To 
illustrate the idea, let us outline the method of attack on the system shown in Fig. 4, 


—=" 
_e 








m6 





Fic. 4. System III. 


the points O, D being fixed, and motion being confined to a vertical plane. This 
system may be regarded as system I interconnected with another system of the same 
type. Thus the primitive may be taken to consist of two systems of the type I, as 
shown in Fig. 5. For each system of type I the metrical tensor and force vector are 
known, being of the form (22) and (23). For brevity, we denote them in shape only 


by the following symbols: 


Pepa f ey [2 } $2). 


Then for the whole primitive system the corresponding quantities are 
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i. ee 0; 0 (1’) 
| ex 2 oe oe 
0 7 0 | @ | | 


which, of course, may be written down at once. The configuration space of the present 
primitive system is the direct product of the configuration spaces of the two systems 
of type I. 


RG 
—s oe 





2a, 








Fic. 5. Primitive system of system III. 


The interconnection of the two systems at B introduces a single constraint, and 
by the method of subspaces the equations of motion of system III may be obtained 
in a routine manner. 

The problems discussed above involved only very simple interconnections. When 
the interconnections are numerous and complicated, and the elements interconnected 
are themselves known complex dynamical, electrodynamical, or hydrodynamical sys- 
tems, Kron’s method of subspaces assumes the highest practical importance. In con- 
cluding the author wishes to thank Mr. Kron for many stimulating discussions of his 
work extending over several years. 

Added April 12, 1944. While the general mathematical theory underlying Kron’s 
method of subspaces as applied to dynamical systems is implied in the “Formal Proof” 
of the present paper, it is not there given in explicit detail. Professor Synge of the 
Ohio State University has suggested that a more explicit proof be included which 
goes directly to the mathematical basis of the method, and has kindly communicated 
the following outline of a method of proof from a different point of view which will 
be of interest to mathematicians wishing to see clearly what is fundamentally in- 
volved mathematically. 

Let (a, b) (e, f) (i, j) take three different ranges of values, with the usual summa- 
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tion convention for each. Let there be two independent holonomic dynamical sys- 


tems: 
I. Coordinates: x¢, 
K. Energy: Ta) =3max%?, 
Generalized forces: Xq. 
II. Coordinates: x°, 
K. Energy: T @) = $mesx*x!, 
Generalized forces: X,. 
Let us define, with D=d/dt, 


Sq = D(OT (1)/0%%) — OT 3) /dx?, 
ea ie ae (A) 
S, = D(OT (2)/dX*) — OT (2)/dx*. 
Then the equations of motion of I and II are Sg=Xa, S,=Xe. 


Now establish constraints between I and II, the reactions of constraint being 


workless. These constraints may be written 


x* = x(x"), x? = x(x), 


where x‘ are the generalized coordinates of the system III resulting from the combina- 


tion. Write 


Cy = 0x°/dx, Ci = 0x°/0 x". 


We have then 


It is easy to prove that 
DCj = ax%%/dx', DC, = d%*/dx’, (B’ 
) 
C? = d%°/ dx", Ci = d%°/ 0x. 


Let Xj, X; be the reactions due to the constraint. We have 


Xiix%+ Xlbx* = 0 


— 
5 
fe) 
= 


for any displacement satisfying the constraints, i.e 


] 

Q 
Q 
4 


bat = Cthxi, bx? 
Hence 
X'C?+ X'Ci = 0. 
Now the equations of motion of I and II under the constraint are 


Sa= Xat+ Xz, S,=X.+X, 
and hence ; 
SaCi + S.Ci = Xai + XG. (C) 
We have for system III: 
Coordinates: x‘, 
K. Energy: Tg) = $mj;x'%', 
Generalized forces: X;. 


Let us define 
S; = D(OT (3) /0x*) —_ OT (3)/0x". 
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Our problem is to show how S;, X; are to be computed in terms of the elements of I 
and II. (We know from dynamical theory that S;=X;, but we can forget this knowl- 
edge, as we prove it incidentally below.) 
We have 
Ts) = Tay + Fe), 


and hence 


ms; = marCiC + myCiC (D) 
It is easily seen by direct transformation that 
S; = SoC + S.C. (E) 
By considerations of work we have 
X bx' = X,bx* + X bx", 
and so 
X;= XaCi + X.C (F) 


Hence by equation (C), we have, as equations of motion of III, S;=X;, where S; and 
X; are given by (E) and (F). The transformation of the metric is given by (D). 

We may sum up essentially by saying: Metric and force transform by (D) and (F) 
when two systems are linked by workless constraints. The extension to the linkage of 
any number of systems is immediate. 
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RIGOROUS SOLUTIONS FOR THE SPANWISE LIFT 
DISTRIBUTION OF A CERTAIN CLASS OF 
ATIRFOILS* 


BY 
OTTO LAPORTE 
University of Michigan 


1. Introduction. The problem of the spanwise lift distribution of an airfoil in a 
uniform stream of air leads, as was shown by Prandtl,' to the following integro- 
differential equation for the circulation ['(y): 


1 ft? dy af ar : 
7} y—y' dy’ sd mc(y)V ™ (1) 
In this equation y is the span coordinate, and —3b)Sy34); V is the velocity of the 
air at infinity, its direction being parallel to the x-axis; c(y) is the chord function de- 
termining the shape or planform of the airfoil; m is a numerical constant, which the 
theory of wings of infinite span fixes at 2x but which seems to have an experimental 
value in the vicinity of 5.5; a is the geometric angle of attack, which for flat wings 
is a constant but for twisted wings is a given function of y. For a derivation, the reader 
is referred to the papers mentioned in footnote 1; we only wish to note here that the 
downwash velocity w is essentially the first term of (1) and has the form 





i 2 +1 dy’ av 


(2) 





w(y) = 


4n 6b J_3 n—7 dy’ 

where 7 = 2y/d is a dimensionless span coordinate. Because of the singular integrand, 
the Cauchy principal value must be taken both in (1) and (2). Of fundamental im- 
portance, furthermore, is the elliptic wing of chord distribution 


c(y) = co/1 — 9? (3) 


for which I'(y) is also of elliptic shape while w(y) is constant. 

In the extensive literature dealing with solutions of (1) for a given planform c(y), 
the following procedure is used almost always: the substitution 7 =sin @ is introduced, 
and the chord function c(y) is developed in a Fourier series. It is assumed that T is 
also of this form, whence equation (1) is satisfied at a discrete number of y values. This 
leads to a system of linear equations for the Fourier coefficients, the solution of which 
is usually extremely laborious. 

In contrast to this method, the point of view adopted by Trefftz? seems to me 
more powerful. This author considers the potential flow in the complex y+iz plane 


* Received Feb. 28, 1944. 
1L. Prandtl, Géttinger Nachrichten, 1918, 451, and 1919, 107. Also N.A.C.A. Report 116, 1921. 


See also the presentation in Mises and Friedrichs’ Fluid dynamics, Brown University mimeographed 


lecture notes, 1941, p. 108 and foll. 
2 E. Trefftz, Zeitschr. f. Angew. Math. und Mechanik, 1, 206 (1921). 
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at a large distance behind the wing. The wake, stretching along the y-axis from 
—b/2 to +6/2 becomes a line of discontinuity along which the velocity potential 
suffers a jump determined by the circulation. The integro-differential equation (1) 
reveals itself as a boundary condition which the velocity potential has to obey along 
the line —b/2 Sy $b/2. “With the aid of conformal transformation the field is brought 
into relation with the field outside of a circle of unit radius; then the potential is ap- 
proximated by a trigonometric expression and an approximate fulfillment of the 
boundary condition is sought.”® 

In the present paper the point of view which was used by Trefftz is extended, but 
instead of “seeking an approximate fulfillment of the boundary condition” for an 
arbitrary chord function c(y), a simultaneous choice of the function and of a con- 
formal mapping function transforms the problem into a boundary value problem of 
classical type, which can be solved rigorously and in every detail. The resulting 
formulae lend themselves readily to numerical computation. 

The family of planforms to which one is led in this fashion is represented by the 
equation 


c(y) = co[(1 — 9%)(1 — x2?) J? (4) 


For 


1, (4’) 


this results in airfoils of taper greater than in the elliptic case (3), while for 


0 


IIA 


K 


IIA 


OS Si, | (4”) 


airfoils blunter than the elliptic ones are obtained. Fig. 1 shows the entire family for 
a fixed span and aspect ratio.‘ 

















Fic. 1: A=8. 


To be sure, for a given planform elaborate approximate methods, such as the 
Fourier series method described above, will probably always have to be used. How- 





3 y. K4rma4n and Burgers, in vol. II of Durand’s Aerodynamic theory, Springer, Berlin 1935, p. 171. 
‘ The aspect ratio cf is defined as b?/S, where S is the area of the wing. For the calculation of S 


as a function of «x, see (22). 
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ever, we shall now show that at least for one planform family, which depends on two 
parameters, the problem can be completely solved. It is hoped that the method used 
here will be generalizable so as to furnish perhaps other chord functions in the future. 
At any rate, for a planform which is close to one of family (4) a method of successive 
approximations can be readily set up. 

It seems that a comparison with the methods used in the theory of wings of infi- 
nite span is not superfluous. Here rigorous solutions are available for certain families 
of simple profiles, the simplest one of which is that of Joukowski. Although the de- 
signer will probably not be satisfied with such a simplification and will turn to more 
elaborate methods, there is an advantage in having easily derived formulae in closed 
form—both as far as quick estimates and classroom presentation are concerned. In 
lectures it is standard procedure to present the flow around the Joukowski and per- 
haps also around the Kaérm4n-Trefftz profiles before turning to more general meth- 
ods; on the other hand, when dealing with wings of finite span, the derivation of 
Prandtl’s equation (1) is usually followed only by a detailed discussion of the elliptic 
case and then by a mere description of the Fourier series methods. The formulae 
given below are intended to fill this gap. 

2. Prandtl’s equation as a boundary condition in the complex plane u=y-+iz. A 
brief recapitulation of Trefftz’s work? is necessary following the presentation of Mises 
and Friedrichs.! In the y+7z plane there exists a flow derivable from a potential 
$(y, z). This potential is everywhere continuous except on the projection of the airfoil 


z= 0, -—-b/2s ys + 3/2, (5) 


where, because of the existence of a circulation, ¢ is discontinuous. Assuming that the 
values of ¢ at opposite points of the slit (5) are equal and opposite, we have 


I'(y) = 4¢(y, 0). (6) 


The downwash velocity w becomes 


dg 
wy) = (~) ; (7) 
Oz 0 
The problem is therefore to find a solution of the Laplace equation 
A¢= 0 (8) 
which on the upper side of the slit (5) satisfies the boundary condition 
dg 8¢(y, 0) 
(“) — —— = — Va(y). (9) 
dz/o = mc(y) 


This condition results from substituting (2), (6) and (7) into (1). The corresponding 
condition on the lower side of the slit differs from (9) in that the sign of the second term 


is opposite. 
Condition (9) characterizes the present problem as a boundary value problem of 
the third kind. As mentioned earlier, Trefftz now maps the region outside slit (5) into 


the interior of the unit circle by means of 


-<(: -) (10) 
 - sla . 








1944] SPANWISE LIFT DISTRIBUTION OF AIRFOILS 235 


Through the introduction of polar coordinates ¢=r(cos 0+ sin @) he is then led to ap- 
proximate solutions having the form of a Fourier series in 0. 

It is evident that, because of the occurrence in (9) of the function c(y) in the co- 
efficient of $(y, 0), exact solutions of the problem will in general be difficult to obtain. 
But let us also examine the effect of a conformal mapping upon (8) and (9). While 
the former is invariant under arbitrary conformal transformations the latter is not, 
because of the occurrence of the first derivative. One may say that each conformal 
transformation causes a new variable coefficient to appear both in the ¢ term on the 
left and on the right side of (9). 

We propose to let these two effects counteract one another and to choose both ap- 
propriate conformal transformations and appropriate chord functions c(y) so as to 
arrive at a boundary condition with constant coefficients which lends itself to treat- 
ment by classical methods. An example will best serve to illustrate this procedure. 

Instead of introducing real polar coordinates within the circle into which the slit 
(5) is mapped by (10), we shall map the interior of the circle onto the strip bounded 
by the points 0, + ©, + 0-+2m7i, 27i of aX plane by means of 


b 
t=e, u = —chi, (11) 
2 
where \=y+7v. The potential must be a solution of the Laplace equation in \ and yp, 


and must satisfy the following boundary condition along the imaginary axis onto 
which the circle | ¢| =1 is mapped: 


Og 45 sin v R 
(*) — ————— (0, v) = — 46Va sin ». (12) 
Ou/o mce(%b cos v) 


In order that the coefficient of @ be a constant we must choose c($d cos v) proportional 
to sin v, whence, by virtue of (11), we have the elliptic chord distribution (3). 

3. A transformation mapping the interior of the circle t| =1 into that of a rec- 
tangle. We now replace the function (11) by others which map the interior of the unit 
circle into various other regions. Each time, in order to have a boundary condition 
with constant coefficients, an appropriate chord function must be chosen. When 
mapping the unit circle of the ¢ plane into a rectangle in the usual Schwartzian way 
such that four arbitrarily chosen points on the periphery correspond to the corners, 
a chord function results which possesses singularities at those points of the span 
which are the images of the corners. It is, however, possible to map a rectangle onto 
the unit circle, such that two opposite sides of the rectangle become two opposite 
semicircular arcs while each of the two other sides of the rectangle map onto a slit 
protruding radially part way into the circle.’ This is illustrated in Fig. 2. The mapping 
function which accomplishes this is 





— | 
t= /k sn —Z, (13) 
Tv 
where the circle is in the ¢-plane (t=r-+is), the rectangle is in the Z-plane (Z=X+/Y), 


’ G. Holzmiiller, Einfiihrung in die Theorie der isogonalen Verwandtschaften, 1882, p. 256 and foll. 
See also Cayley, Collected papers, vol. 13, papers No. 891, 920, and 921. 
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k is the modulus of the elliptic function and is between 0 and 1, K is the complete 
elliptic integral. We now separate real and imaginary parts by means of the addition 
theorem of the sn function and put 








ee (14) 
“£2 
In view of the fact that sn(tK’/2) =zk-'/?, we obtain readily 
sn 2KX/x cn 2KX/xdn 2KX/xr 
r= (1+ b) s= + (15) 
1+ ksn? 2KX/x 1+ ksn*? 2KX/xr 
Therefore 
r? + 5? = 1, (16) 


Thus the upper straight line (14) maps into the upper semicircle, and lower one into 





























t- plane z-plane 
s y 
fe (i) 
a 
aK 
4K 
fe i) fe | i. 
Kk] a 
2 
(4} 
Fic. 2. 


the lower. If in (15) we set X = + $7, we obtain s =0. The four corners of the rectangle 
are hereby fixed. To see what corresponds to the vertical sides, we set Z= +7/2+i1Y 
in (13), to obtain 
2K 
r=+v¥knd(—Y, #’), (17) 
T 
where k’?=1—k?. The slits protruding into the circle along the real axis therefore 


terminate at r= +k"/?, 
4. The boundary condition in the Z-plane. For the upper side of the rectangle this 


condition becomes, by use of (9), (10) and (13), 





0g 8 du a K’ du 
=) iene =) (x, = ) . =) (18) 
oY 4K’ /4K mc dZ 4K’ /4K 4 K dZ eK’ /4K 
where 
du b a 2K 2K 2K 
= VEK(1- ns — Z) en——Z dn—Z (19) 
d 2a T T T 


and Z must be put equal to X+irK’/4K. The chord c(y), now regarded as a func- 
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tion of X, must be chosen proportional to (19). Therefore it becomes necessary to 
express (19) as a function of y in the original u plane. We have at first from (13), 


dus b/RK 1 i? 1/2 
wo (AC e-w]" 
dZ Qn i? k 


where | z| = 1. Using for the moment the angle y as in (12), this may be written in the 


form 

du ib/ kK 1\7/ 

—_ = sin | 20082 - (4+ —)| 
dZ T k 








which, since cos y=, becomes finally 


du 1 
— = — bK(1+ &)[(1 — n2)(1 — wm’) ]2. (20) 
dZ T 
The abbreviation 
~= (+ vi) (21) 
kK = 2 Vk 


was used, where x is between 0 and 1. We see thus that we are led to the planforms 
(4), (4’) which represent wings of taper greater than the elliptic wing. 

Before (4) and (20) are substituted into (18), it is convenient to calculate the area 
S and the aspect ratio -4 =6?/S. It stands to reason that the area 


6/2 
S = f co[(1 — ?2)(1 — xn?) ]"2dy 

—b/2 
is expressible in terms of the complete elliptic integrals E and K of x, but in view of 
the fact that (21) is the well known Landen transformation, S can also be reduced to 
the complete elliptic integrals of the modulus k. The result is 





S = beG(k), (22) 
where 
1 i+ 6k+ 
a = [7 aw - a - ba - 30x]. 22 
(®) = |g B® ~ A= Od = 3K (221) 
G(k) is a purely numerical constant depending on taper. For the aspect ratio we have 
. b 
= , (23) 
G(k)co 
and for the average chord é, 
¢= G(R)co. (23) 


Substitution of (4), (20) and (23) into (18) gives the final form of the boundary 


condition, 
od 8 A a K’' 
() +— —(1+ WGKe (x, = =) 
K’/4K 4 K 


oY x m 
ia 2K 2K 2K 
=—y BK [a(1— &*nst —Z)en—z dn —zZ . (24) 
2x 7 T w YerK’ /4K 
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Variable coefficients now only occur on the right. It should be noted that for 
Z=X+inK’'/4K the right-hand side represents a real function of X. 

5. Solution of the boundary value problem for tapered wings. This solution has 
to satisfy the following requirements: 

(1) Inside the rectangle |X| $7/2, | Y 
tion. 

(2) On the top side Y=7K’'/4K it must satisfy (24). 

(3) On the lower side it must satisfy a similar condition differing from (24) only 
in the sign of the first term. 

We shall seek a solution of the form 


<7K’'/4K, it must satisfy Laplace’s equa- 





@(X, Y) = pS sh nY(A, cos nX + B, sin nX). (25) 


n 


In this series each term is a solution of the Laplace equation. The sh function of Y 
is alone present because of condition (3), or in other words, because we want ¢ to 
have equal and opposite values at opposite points on slit (5) in the y+7z plane. To 
determine A, and B,, the right-hand side of (24), considered as a function of X, has 
to be expanded in a Fourier series. 

At this moment it should be remembered that for twisted airfoils a is a function 
of y or 7. It will be an even function as long as the ailerons are in their normal position. 
Otherwise, a may have an odd component or may even be discontinuous. At the pres- 
ent time only the case of constant a@ will be treated.® 

To develop the right-hand side of (24) into a Fourier series we begin with the 
familiar series’ for sn: 


2K = : T % 
sn —Z = be a, sin (2n+1)Z, a, = a cosech x(n + 3) 3 . (26) 


T n=0 


This series converges uniformly as long as | Z(z)| <7K’'/2K, i.e., within a horizontal 

strip whose median line is the real axis and whose width is rK’/K. But since this is 

just twice the vertical dimension of the rectangle of the boundary value problem (see 

Fig. 2), series (26) will certainly converge absolutely anywhere that it is needed at 
present. We differentiate and put Z=irK’/4K +X, to obtain 

2K 2K . = _« K’ : 
en —Zdn—Z = — )) (2n + 1)a, cos (2 4+- 1)( i— —+X). , 
c 4 K (27) 


T T 2 n=0 





To get the other part of the right-hand side of (24), we put Z=irK’'/2K+2Z’, 


whence 
2K 2K 
ns —Z = ksn—Z’. 
T T 


6 More general cases of twisted airfoils are being computed at present and will be reported on else- 


where. 
7 See, for instance, Whittaker and Watson, Modern analysis, 4th ed., Cambridge University Press, 


p. 510. 
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The application of (26) with Z’ as the variable gives 


2K ~ 
k-! ns —Z = D> ay sin (2n + 1)Z’. 


T n=0 
After differentiating with respect to Z or Z’, we put Z’ = —i7K’/4K+X (which is 
the same as putting Z=+irK’/4K+X), to obtain 
2K 2K ri 
- kia—Z2 dé —éioe-— 


T T n=0 








- (2n + 1)a, cos (2n + »(- i - + x). (28) 


The sum of (27) and (28) furnishes the desired series which, because of the second 
formula of (26), may be written in the form 


cos (2n + 1)X 


2K 2K 2K a 
(: — k™ ns? — z)en—Z dn— = - >> (2n + 1) » (29) 
T T tg 2K7R amo sh (n + 4)xK’'/2K 








where it is to be understood that Z =irK’/4K+4X. This series converges uniformly 
for all real values of X. 

The final solution of the problem can now be written down. Introducing (29) and 
(25) into (24) one obtains for the velocity potential, 











Vab 2 sh (2n + 1)Y cos (2n + 1)X 
$(X, ¥) = — —_}> ri he al (30) 
4 KVk Dash (2n+ 1)xK'/4K 
where 
D, = ch (2 4oe. 0g ISOS a pans (31) 
.= 2 — —+— — — — sh (2n — . 
gas Sr are ee se 4K 


This series converges uniformly as long as | Y| <mK'/2K, which is more: than we 
need for numerical calculations, since the rectangle only extends to Y= +7K’'/4K. 
Of greater practical importance than ¢ is the sectional lift or the sectional lift 
coefficient c; which may be defined as follows: let / be the lift per unit span, 


L = 4pV%éc1, (32) 


where é is the average chord (23). Since from (6) / is also 49 V¢(X, 7K’/4K), the 


following expression results: 





= 4 
C1 = 2ra _ be — cos (2m + 1)X. (33) 
Kv k n=0 n 


For this rapidly converging series, numerical tests have shown that three or four 
terms suffice to give a result correct to one part in ten thousand. To formula (33) 
should be added the relation between X and the original span coordinate y. From (15) 
and (10) it is found to be 

sn 2KX/x 
1+ ksn? 2KX/xr 





y= “ (1+ k) (34) 


Formulae (33) and (34) give the spanwise lift distribution as a function of y. 
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6. Results for the blunt wing. The method of the previous section may readily 
be used to get similar results for wings blunter than elliptic ones. The underlying idea 
is the following: in section 2 we mapped the ¢ circle onto a square such that the semi- 
circle in the upper half plane became the upper side of a rectangle in the Z plane, 
and correspondingly for the lower semicircle. Now, the interior of the unit circle in 
the ¢ plane will be mapped so that the right and left semicircles will correspond to the 
right and left sides of the rectangle, while two slits from +i to +i+/k will correspond 
to the upper and lower sides. In short, the role on the real and imaginary axes in both 
planes will be reversed. Results will merely be given and formulas will receive the 
same numbers as the corresponding formula of the “tapered” case, except that primes 
will be added. 
The mapping just described is performed by 


> 


t = i-4/k sn — iZ. (13’) 
T 


Upon putting X = +7K’/4K the point in the ¢-plane moves on a circle r?+s*?=1, 


























t-PLANE Z-PLANE 
Y 
S 
[3] (2) [1] 
(3) } (1) 
 s 
2 
(2) 
r x 
ts) Le od 
4K 
(4) | (6) 
[4] (5) (6) 
Fic. 3. 


while Y=+72/2 gives the slits along the imaginary axis of ¢. This is shown in Fig. 
3. The interior of the unit circle of the ¢ plane is then mapped by means of (10) onto 


the exterior of the slit (5) in the u plane. 
The boundary condition, now along the two vertical sides, becomes 


a¢ 8 /du x K’ du 
=) +— =) o¢{— —> v) = 11Va =) : (18’) 
OX/ exyax me \ dZ / ex sk 4K dZ / ek’ /4K 


with 
du b 2K 2K 2K 
=) = —VEK(1 + k-' ns? — iz) cn —iZ dn — 72Z, (19’) 
dZ / ek! 14K 


2a T rg T 


which is imaginary for Z=2rK’/4K+iY. The determination of the planform results 


from expressing this in terms of 7: 
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d 
a(S) = eK - HVE ITE, (20) 
dZ / ek’ /4K T 


where 


L-H(G-u) ev 


The different signs should be noted. The two ;o9 
functions x of x’ of k are plotted in Fig. 4. 
In this way a planform blunter than the elliptic 
one results, namely (4), (4’’). When x’ becomes 
greater than one, the chord function has a maxi- 
mum between 7=0 and 1. The bluntest wing .6 
shape without such a maximum is attained for 

x’=1 or 4 


kor =3- J/8, (35) 





8 





in which case the lemniscate functions result. 
Thus the formulae of this section, while being 
valid for any k between 0 and 1 have aero- O 
dynamical importance only for k between 0 
and kcr.® 

The area of the wing becomes 











S = beoG(— k), (22’) 


where G(—k) is the function in equation (22;) for negative k.® In the case of (34), 
which interests us especially, the elliptic integral becomes reducible to Gamma func- 
tions: 


G(— hen) = TG/2)1(5/4) = .87402 (36) 
J CR T(7/4) ° be 


Formulae for -4 and the mean chord é are taken over with G(—). 

For a constant angle of attack, the only case treated below, the function in (19’) 
will be developed in a Fourier series on the vertical sides X = +7K’/4K of the rec- 
tangle. This is readily accomplished, the difference from (29) being that only terms 
in sin (2n+1) Y occur. The final result for the velocity potential is 


x Vab 2 ch(2n+1)X sin (2n + 1)Y 








¢(X, Y) =— . (30’) 
4 KVk > A, ch (2n + 1)rK’/4K 
where 
« K’ 8A(1i— &#G(— AK x K’ , 
Bb. = tc (In + 2) 4 — — ch (2n + 1) — —- (31’) 
) 4K «fm n+1 ( 4K 


* The mapping for k =3—4/8 is studied and illustrated by figures in Holzmiiller’s book (see foot- 


note 5). 
® However, E and K depend only on ?. 
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Convergence properties are the same as before. The sectional lift coefficient (32) is 
given by 

= > ~ si 2n + 1)Y¥ (33’) 
ee —- gin (2 ; 

Ky k n=(0 A, 





ci = 2ra 


while the relation between the variable Y and span coordinate y is now 


b cn2KY dn 2KY/zr 
_.—— ‘ (34’) 
2 1+ ksn? 2KY/r 





7. The elliptic wing as a limiting case. From (20) or (21) it is evident that the 
chord fur -tion represents an ellipse when k goes to zero. For decreasing k the elliptic 
integrals K and G approach 7/2 and 7/4 respectively, while 


4 
lim (x — log =) = 0, 
k0 k 


The quotient K’/K becomes large and the convergence of the series (33) improves. 
Keeping only the first term we see that 


; ch\ x K’ 
m ()EE-r}-e 
k-0 sh/ 4 K 


4 cos X 
ci = — am 


T 1+ xeA/m 


Thus 
(37) 


which is in complete agreement with Durand II, p. 169, since now cos X = (1 —7n?)*/?. 

8. The limiting case of extreme taper, k =1. This case is of more interest than the 
preceding one, inasmuch as it leads to new and simple formulae in closed form. It 
follows from (20) that the planform is now given by 


¢ = ¢o(1 — 7), (38) 


which represents a parabolic arc or two such arcs joined along the y-axis (see Fig. 1). 
Although the lift distribution can be obtained from (33) by letting k tend towards 1, 
we prefer to derive it afresh. 


The transformation 
t= tanhZ (39) 


maps the interior of the unit circle in the ¢ plane onto a strip parallel to the X-axis 
and bounded by Y= +7/4. The points t= +1, where the ¢ circle crosses the real 
axis move, respectively, towards + © and — ©. Transformations (10) and (39) to- 
gether map the outside of slit (5) in the u plane on the above mentioned strip, in such 
a way that the upper and lower sides of the slit become, respectively, the upper and 
lower sides of the strip, while the tips (y= +5/2) move to infinity. 

The boundary condition on the upper side is the same as equation (18), except 
that the subscript on the various derivatives now is Y =7/4. Since the transformation 
from the u to the Z plane is now 


u = $b coth 2Z, (40) 


as is seen from (39) and (10), we have as the analogue of (19) the equation 
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du 
(=) = b sech? 2X, (41) 
dZ Yur/4 
which, when expressed as function of 7, becomes b(1—7?). Hence the chord is chosen 
as in (38). We note the formulae for mean chord é, area S, and aspect ratio R: 

2 2 3 5b 

€ = — 6; S=—bo; A=——- (42) 
3 3 2 


Co 


The final form of the boundary condition on the upper side of the infinite strip is 


0d 16 4 T 
(=) Pe. “o(x, =) = Vab sech? 2X. (43) 
Vien 3m 4 


Due to the fact that the range of X is from — ~ to + ©, the right-hand side of (43) 
must now be expanded in a Fourier integral. This can be done by standard methods, 
the result being 


° td x 
sech? 2X = — Len Beet. Fava (44) 
0 sh ~ sh xt/4 


For the velocity potential ¢(X, Y) we assume a solution of the form 
¢(X, Y) = f dtA(f) sh ¢Y cos ¢X, (45) 
0 


where the integrand satisfies the Laplace equation. After introducing (45) and (44) 
into (43), A(€) is easily determined. The final form for ¢(X, Y) is 


’ 1 * shi? cota 
¢(X, Y) = — vab f — df, (46) 
4 o D(s) sh xf/4 


where 
16 cA sh f/4 
D(g) = tnt $= = See (47) 
3 m t 


The sectional lift coefficient assumes the form 
= cos ({X 
c1 = 2cfa f dt. (48) 
0 Ds) 


The integral in (46) converges for all X and Y within the strip under consideration. 
In order to obtain from (48) the lift distribution as a function of y, we note here that, 
due to (40), 





y = 36 tanh 2X. (49) 


The above integral may be brought into a rather different and interesting form 
which makes the numerical evaluation very much easier. Let us consider the integral 


| i rdr cos té (50) 
aie chester. 





which, except for a simpler notation, is the same as that in (48). On the imaginary r 
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axis, the integrand has infinitely many poles, which can be computed as the roots of 
the transcendental equation 


ocos¢+asine = 0. (51) 


Denoting the mth root by @n, we find from this that 
cos o, = (— 1)"a(a, + o%)—/2, sin ¢, = (— 1)"*e,(a? + 02)-"/2, (51’) 


For r=0 the integrand is regular. It now becomes convenient to decompose the cos 
in the numerator into two exponentials. Two integrals arise in this way, one contain- 
ing exp irt, the other exp (—iré). According to a standard argument, the integral with 
the positive exponential tends to zero when integrated along the circumference of a 
large semicircle in the upper half plane, as the radius increases; similarly for the other 
integral and a semicircle around the origin in the lower half plane. Thus the integral 
with the positive exponential, taken along the real axis from — © to + ~, is equal to 
the sum of the residues at the poles ¢, which lie in the upper half plane, and the in- 
tegral with the negative exponentials is equal to the sum of the residues at the poles 
of the lower half plane. 

Applying these considerations to (48) we see that the quantity a, which deter- 
mines the roots of the transcendental equation, has the value 


a= — — (51”) 


while the sectional lift coefficient is given by 


2 o,[a? + on, |*/2 
ci = Bac >> (-— 1)**! 
n=O a(a+ 1) + 03 
This is a very rapidly converging series, which, although especially convenient near 
the wing tips, may be used to advantage for values of 7 as small as .3 or .2. Very near 
the wing tips X is approaching ~. Hence, using (49) we may write approximately 


a(a+1)+ 0; 2 


which shows that c; has a vertical tangent at the tips since o; is between 7/2 and 7. 
The only inconvenient region is that of small X, especially the point X =0. It 
seems that only graphic integration can be used to find c; at the center of the wing 
from formula (48). 
9. Total lift and total induced drag. For these important quantities rapidly con- 
vergent expressions are readily derived. The total lift L is obtained by integrating 


pVT across the span, 


e~tnX/t, (52) 








c= 


+b/2 
ica f $(y, Ody, 


6/2 


and transforming to the Z plane, to obtain 


+2/2 Tr P du 
L = 4pV 6(x.= —)\(S) dX. 
—s/2 4 K dZ / ek" 4K 


2 
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The two Fourier series for @ and du/dZ are then substituted from (30), (19) and (29). 
From the completeness relation it follows that the integral of the product of the two 
series is equal to the sum of the products of the coefficients. Instead of the for- 
mula for ZL, merely the final formula for the over-all lift coefficient cz, defined by 
L=$pV*Sc_z, will be given :!° 

nr Aa & 2n+ 1 

4 K*k j20 D, sh (2n + 1)xK'/4K 


The abbreviation D, was defined in (31). In the limit of vanishing k this reduces ex- 
actly to formula (2.15), p. 169 of Durand, Vol. IT. 
For wings blunter than elliptic ones the result is 


‘ds 2n + 1 
65 tae > as “ 


4 K*k < A, ch (2n+ 1)4K'/4K 


where A, is the coefficient defined in (31’). 
For the parabolic wing treated in section 7, it is necessary to go through the opera- 
tion just described for the Fourier integrals (44) and (46). The lift coefficient is now 


T ces d 
gale f a. Sn (54) 
Z o D(g) sh rf/4 


For the induced drag D; similarly simple formulae may be obtained. Using (7) 
we have first in the original u-plane (u=y+7z) 


(53) 





(, = 





(53’) 





6/2 0g 
D; = 4 , 0 — d . 
p o(9 (=) y 


—b/2 


In the Z plane of the tapered wing this may be written 


of x K’ 0p 
D; = 4p f of 25 =) (—) aX. 
—s/2 4 K OY / eK’ 4k 


Now the completeness relation must be applied to the Fourier series (30) and its 
own derivative. The result is, expressed in terms of the drag coefficient, 











wm arA © In+1 h (2m +1) a K’ (55) 
=—_— cot n ve ’ 
° 4 Ke Dt 4K 
while for the blunt wing we have 
* ak < 2nt+1 x K’ 
Cp = wi ps tanh (2m + 1) — —-» (55’) 








_ 4 K*k,2o A? 4K 
10 Formulae (53) and (54) may be derived in a different way without appealing to the completeness 


relation. The velocity potential ¢(y, s) is only the real part of a complex stream function F(u) in terms 
of which the total lift may be written in the form 


L = 2pVRe $ F(u)du. 


When transformed into the Z-plane, this becomes a contour integral around the rectangle of Fig. 2. 
This integral can be carried out simply by applying Cauchy’s theorem. The same can be done for the 
rolling moment. This method recalls the Blasius theorem in the theory of infinite span. 
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and for the parabolic wing 


™ ° $dg ug 
cp = at f —— coth — (56) 
2 0 D*¢) 4 
Formulae (53) to (56) give the dependence of lift and induced drag upon aspect 
ratio c4 and taper k. However, it is not necessary to compute more than c_, for there 
exists the following simple relation between c, and cp: 





OCL 
cp = acA ; (57) 
dcA 
Thus after having plotted a set of curves showing the dependence of cz, upon the as- 


pect ratio, one may obtain cp by graphic differentiation. 
Further over-all quantities of practical interest are the rolling moment 





b/2 
tov f o(y, O)ydy, 


—b/2 


4 fo » (=) d 
$(y, 0) (—) yay, 
4 — OZ eg , 


but since in this paper we have confined ourselves through out to constant angles of 
attack, they will vanish. However, when a(y) is assumed to have an odd component 
the methods described here still work and the above moments yield formulae of type 
(53) to (56). 

10. Numerical results. All formulae given in the preceding sections for the sec- 
tional lift coefficients are readily evaluated, due to the fact that the series converge 
rapidly. An exception is c; for the parabolic wing shape for small or zero span coordi- 
nate, in which case it is necessary to resort to graphic integration. Following 
v. K4rm4n and Burgers,’ we plot c;/ma as a function of the span coordinate 7 rather 
than c¢; itself. In this way it is necessary to vary only the two parameters¢e/4/m and 


and the yawing moment 


the taper constant k. 
In Tables 1 and 2 values of c;/ma are given which are computed from formulae 


(33) and (34) for k=+/.1 and k=+/.2. In the first column are the values of 2KX/z 
which quantity serves as a convenient parameter, in the second X, in the third 7, 
computed by means of the tables by Milne-Thomson". In the last three columns the 
values of c;/ma are to be found for -4/m equal to 1.0, 1.5 and 2.0, corresponding to 
aspect ratios of about 5.5, 9.25 and 11. 

Among wings blunter than elliptic only the one with keg of (35) has been com- 
puted using formulae (33’) and (34’). However, since elliptic functions with k values 
as small as kon were not tabulated by Milne-Thomson, equation (34’) was expanded 
in powers of that small constant. Values of c;/ma are given in Table 3. 

Finally in Table 4 the same data are given for the parabolic wing of section 7. 
In column 1 are the values of the parameter X, in column 2 the 7 values obtained from 
it with (49), and in the next three columns the values of ¢c;/ma. For X =0, .1, .2 they 
were obtained by graphical integration from (48), and for all greater values of X from 


11 L, M. Milne-Thomson, Die elliptischen Funktionen von Jacobi, Springer, Berlin, 1931. 
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(52). The roots of the transcendental equation (51) do not seem to have been calcu- 


lated before, at least not for the values (51’’) of the constant a. They were determined 


to five significant figures. 
These tables are supplemented by the following figures: Fig. 5 makes possible a 


A 


l4r 


k={02 


= 
" 
oe 
Nw 








Fic. 5. 


quick determination of the span coordinate associated with a certain coordinate value 
on the side of the rectangle of the boundary value problem. The two curves marked 
k=4+/.1 and k=+/.2 represent plots of equation (34) and go with Tables 1 and 2. 


TABLE I 








Tapered Wing k=+/.1 











| C1/me 

2KX/x x | n | cA/m=1 cA/m=1.5 | cA /m=2 
0 0 | 0 1.0191 1.1207 1.1807 
1 09742 13096 1.0070 1.1068 1.1655 
a 29225 37801 91652 1.0026 | 1.0527 

5 .48709 .58643 76447 82908 | 86553 

8 77934 | 80778 +| ~~ $0527 53846 55574 

1 1.07159 |  .93247 | —.28109 29429 | 30021 








The third curve is a plot of equation (34’) and goes with Table 3. Figs. 6, 7, 8 show the 
dependence of lift coefficient upon aspect ratio and taper. They should be compared 
with Fig. 75 of Durand, vol. II (which was computed by the approximate method de- 
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Fic. 6. A /m = 1.0. 














Fic. 7. Al /m = 1,5. 
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Fic. 8. c4/m=2.0. 
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TABLE II 
Tapered Wing k= 4+/.2 





















































| ci/ma 
a — — es | 
2KX/x | x | » | A/m=1 | cA/m=1.5 A/m=2 
ee een ee | 
0 0 | 0 1.0306 | 1.1365 1.1987 
1 09465 | 14380 1.0150 | 1.1184 1.1792 
3 28394 | 41052 90134 | 98667 1.0370 
5 47324 | 62545 .72028 .77905 .81219 
8 75718 | 83613 44105 46565 .47700 
1.1 1.04113 | 94328 .23139 23804 24100 
TABLE III 
Blunt Wing k=3—1/8 
C:/ma 
Y ” A/m=1 cA/m=1.5 A/m=2 
x/2 | 0 .91288 .98248 1.0205 
1.2 .30920 88828 96045 1.0006 
1.0 47438 84591 91985 96191 
7 70700 71481 .78512 82654 
5 83962 | 56705 62678 66262 
35 91816 | 42178 .46800 .49601 
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TABLE IV 
Parabolic Wing k =1.0 

















c1/ma 
X n 0A /m= cA /m=1.5 | 0A /m=2 
0 0 1.050 1.164 | 1.230 
i | .19737 1.020 | 1.130 | 1.166 
2 .37995 .9297 1.029 1.057 
4 66404 66568 .71789 . 74320 
6 .83365 .42208 .43141 .44003 
9 .94681 .17295 .17296 | 17297 


scribed in the introduction). The lift coefficient shows, of course, the expected span 
dependence, inasmuch as a more highly tapered wing gives rise also to a more highly 
tapered lift curve. Of special interest is the parabolic case, for which the lift coefficient 
goes to zero with infinite tangent although the chord function c(y) has a finite tangent 
at the tips. 

The author should like to thank Mr. H. Yoshihara for extensive assistance with 
the numerical calculations. A grant from the Faculty Research Fund of the Univer- 
sity of Michigan is gratefully acknowledged. 
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SOLUTION BY RELAXATION METHODS OF PLANE 
POTENTIAL PROBLEMS WITH MIXED 
BOUNDARY CONDITIONS* 


BY 


L. FOX 
Imperial College of Science and Technology, London 


1. Introduction. The method of relaxation, as originally propounded by South- 
well [1], was used to calculate the stresses in braced frameworks. A physical picture 
of the method, as presented by him in the Wright Brothers Memorial Lecture for 
1941, is the following. At each joint of the structure constraints are applied which 
prevent joint displacements and bear all the load. One constraint is then relaxed, 
thereby transferring some of its load to the members of the framework and some to 
adjacent constraints. Each constraint is relaxed in turn, and more of the load is 
imposed on the framework, until the residual loads (still borne by the constraints) 
may be deemed negligible. 

In a series of eight papers [2:I-VIII], Southwell and collaborators have applied 
relaxation methods to various engineering problems. In some of these the method is 
applied to two-dimensional problems [2:III], and solutions are obtained of the 
equation 
Vw=Z (1) 
for any boundary on which w is prescribed, Z being a given function of x and y. 
Here Prandtl’s membrane analogy [3] is used, in which w is the displacement of a 
membrane fastened at its boundary and acted upon by a transverse force Z. The 
membrane is replaced by a mesh of uniformly tensioned strings, the mesh lines 
forming squares or equilateral triangles, and the tension in the strings being propor- 
tional to the surface tension of the membrane. Initially the mesh is flat and the 
load Z is taken by constraints acting at the mesh points. The constraints are relaxed 
one by one, just as in the framework, until the loads are all taken by the strings, and 
the resulting displacements of the mesh points are recorded. Evidently, as the mesh- 
length decreases, the approximation of the mesh to the continuous membrane is im- 
proved. 

In a recent paper [2:VIII], Southwell and Vaisey extend the membrane-net 
analogy to obtain solutions of Laplace’s equation in the case when the normal 
gradient dw/dv, instead of the function w, is given at the boundary. Here 0w/dr is 
regarded as a line intensity of transverse loading applied round the boundary of the 
membrane. This load is then integrated and distributed statically to the strings 
which cross the boundary. 

There is a mathematical treatment of the above problems, based on finite differ- 
ences. For Laplace’s problem of the first kind, in which the function is specified on 
the boundary, the finite difference equations have been derived [2: III]. They are in 
general identical with the equations obtained from the analogy of membranes and 


tensioned nets. 


* Received May 19, 1944. 
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The corresponding mathematical treatment of Laplace’s problem of the second 
kind, however, in which the normal gradient of the function is given on the boundary, 
was not given by Vaisey and Southwell, and it is this treatment with which this 
paper is concerned. The finite-difference equations are substanfially different from 
those obtained in their paper. 

The technigue presented here would seem to be particularly desirable for prob- 
lems in which the boundary conditions involve both the value of the function and its 
normal gradient. For then the mechanical analogy becomes somewhat complicated, 
especially in the case of solids of revolution, for which the analogy of variably ten- 
sioned nets is not attractive. 

2. The finite difference approach to the relaxation method. For square meshes 
(Fig. 1), we have the approximations 





aw , ow 
2a —— = WwW — Ws, a*—— = wi + ws; — 2w0, 
Ox Ox? 
a’v?w = Ww - We + W3 a Ws 4wo. (2) 


Similar formulae are easily obtainable for triangular meshes. The order of the error 
in these approximations is known in each case, and decreases with the mesh-length a. 
There is a finite-difference equation of type (2) for every mesh point, and the 
solution of these equations gives a numerical value of w at each mesh point. 
At points close to a curved boundary, such as 0 in Fig. 2, we obtain the finite- 
difference equation as follows, in the case when w is given on the boundary. The 






































‘ 
2\ x 2 
—_ 
3 O B | Cc F/ O " 
ane) 
Fic. 1. Fic. 2. Fic 3. 


arm 01 passes outside the boundary, and a linear interpolation along 01 yields 
Wp=Wp+(w1—Wo)h, from which we obtain for the finite-difference equation at 0 


We + Ws + wy + —— (s + ~) wo = 0. (3) 
h h 

All these formulae are reproduced identically by the net analogy, but a more 
accurate formula than (3) given by Christopherson (4), and obtained by a parabolic 
interpolation, has not been deduced by analogy. 

When 0w/dr, rather than w, is given on the boundary, the same general equations 
hold as before, but there is a different procedure for points adjacent to the boundary. 
To write down the finite-difference equation at the point 0 (Fig. 3), we require 
wa, and wg. For wa, we draw the normal APE and use the approximation 
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Ow 
AE: (=) = Wa — We. (4) 
Ov Jp 


A linear interplation along OG gives 
OG: we = OE: we aa EG: wo, 
and elimination of wg between these two equations yields 


OE 4 +aE(—) 
Wa = — we + — —). 
aa a av |p 


Similarly, we obtain wg by drawing the normal BOF. 

The normal can be terminated on any convenient line. Thus in Fig. 3 we could 
produce AE to K on OD, to obtain wy, in terms of wo, wp, instead of wo, we. The 
shorter the normal, however, the more accurate is the 
approximation (4), so in this case E is the best place 
to stop the normal. On the other hand, the normal 
from B is continued to F, because a termination on 
the diagonal HO would involve the value wz, itself > 
fictitious, and the calculation would become more c 


y 











cumbersome. 

The approximations employed to date have as- \ 
sumed w to be linear along any line near the bound- 
ary. Improvements in accuracy can be made at the 
cost of additional labour. Thus, if the normal is 
stopped so that it is bisected by the boundary (Fig. 4), Fic. 4. 


the formula 

Ow 

AC: (=) = wa — We (S) 
Ov Jp 





assumes a parabolic variation of w. The point C in general no longer lies on a diagonal 
or mesh line, but its value can easily be found, by double interpolation, to the approxi- 
mation of Eq. (5). This procedure yields a higher accuracy, but it is better, except 
when a very high accuracy is required, to use a linear variation together with a finer 


mesh. 
3. Problem I. Let us find the function w, harmonic in the circle 


x?+ y? — 2x — 2y+1=0, 
and satisfying the boundary condition 


df. (6) 


Ov x? + y? 


This is one of the problems attacked by Vaisey and Southwell [2: VIII]. It has the 


exact solution 
w = tan“ y/x, 
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whence we have a guide to the accuracy of our method. The solution is unique, except 
for an arbitrary constant, and we choose the constant so that w=7/4 at the centre 
of the circle. 

First, we take the mesh-length equal to the radius of the circle. The mesh con- 
tains only five points. Fig. 5 shows this mesh, the finite-difference equations, and the 
values of w multiplied by 1000. External mesh points are denoted by open circles. 
For comparison, the exact values are entered under the approximate values. 

When the mesh-length is halved, there are twelve fictitious mesh points and 
thirteen simultaneous equations. The solution is given in Fig. 6. 














7 
1096 @) 
107 
2(5) + .5(2) + .5(4) — 3(1) — 733 =0 (1) 
\ 2(5) + .5(1) + .5(3) — 3(2) + 733 = (2) 
_— = = te 9 ea @ od 2(5) + .5(2) + .5(4) — 3(3) + 2333 = 0 (3) 
> 2(5) + .5(1) + .5(4) — 3(4) — 2333 = (4) 
(1) + (2)+ (3) + (4) —4(5) = (5) 
e } 
0 
Fic. 5. 


As the form of dw/dy indicates, w is skew-symmetrical about the line x =y. This 


feature was not utilized in construction of Figs. 5 and 6, but it is found 
labour-saving device in the case of a mesh-length of one-quarter of the radius. The 


useful as a 


if 
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Fic. 6. 
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results in this case are shown in Fig. 7, and are very close to those for the exact 
solution, the error being greatest at points nearest the origin. This was to be expected, 
since 0w/dv changes rapidly in value across the line x = y when x*?+? is small. Com- 


parable errors were found in the treatment by Vaisey and Southwell. 
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An interesting point is the oscillatory nature of the convergence of values of w. 
In the boundary-value-specified problem, values usually converge from one side 


only. 
4. Problem II. Let us find the function w, harmonic in the same circle as before, 


but satisfying the boundary condition 


0 1 1 y 
(—-—)w-5 y-x-rtant2), 
Ov r r? x 


where r=+/x2+y*. The exact solution is again w=tan y/x. 
Here the boundary condition involves the value as well as the normal slope of 
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the function, but only a slight extension of the method of the previous problem is 
needed. As before, fictitious points are eliminated from the governing equation by 
means of the boundary condition. 

In order to write down the finite-difference equation 
at 0 (Fig. 8), we need w, and wz. As before, the normal AF 


yields 
ow 
AF.- (~) = WA — Wr, 
Ov Q 


and linear interpolation on OE and AF gives 
OF -wr = OF -wz + FE- wo, 
Fic. 8. AF- we = AP: wr + PF: wa. 





From these three equations and the boundary condition (7), wa can be found in 
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terms of wo and wg. A similar operation yields wz, and hence the finite-difference 
equation for the point O can be written down, and the problem solved. 

As in the previous problem it was found that decrease of the mesh-length results 
in oscillatory convergence. The oscillation is rather more violent, but the final result 
shown in Fig. 9 is no less accurate, notwithstanding the additional interpolation. 

5. Summary. In this paper solutions are obtained of Laplace’s equation with 
boundary condition involving either the normal gradient only, or both the boundary 
value and the normal gradient. The need for a technique for problems of this kind 
has arisen in recent developments of the relaxation method. Two problems are 
solved, both have a known analytical solution, and good results are obtained in 
each case. The method used is independent of the analogy of tensioned nets, and can 
be applied without modification to problems for which analogies may be difficult to 
use. 
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—NOTES— 


THE METHOD OF STEEPEST DESCENT FOR NON-LINEAR 
MINIMIZATION PROBLEMS* 


By HASKELL B. CURRY (Frankford Arsenal) 


1. Introduction. The problem considered here is that of minimizing a function of n 
real variables, G(x1, - - - , X,). The object is to find a practical method for evaluating, 
approximately at least, a stationary point for G. 

This problem includes as a special case that of solving a set of simultaneous equa- 


tions 
fil%i1,+++, Xn.) = 0 (4 = i; 2,:-+,m), (1) 


because the function 


G(a1,°-+,%) = ft (2) 


k= 


has a minimum at a solution of (1). It also includes that of determining the parame- 
ters x1,---, x, of a function f(u; x1, - - - , X,) so as to get the best approximation, 
in a least square sense, to a function F(u) for certain values of u; the G in this case 
is of the form given by 


Pp 
G(41,°°*,%&) = > [F (ux) — f(r; %1,°°°, Xp) |. (3) 


k=1 


Certain engineering applications of the latter sort of problem arose in the work 
of the Engineering Research Section, Fire Control Design Division, at Frankford 
Arsenal. In these applications, the function f(u; x1, - - - , x») was sufficiently compli- 
cated so that the standard method for dealing with non linear least square problems! 
failed to converge. Two techniques for dealing with this situation were developed by 
the section under the direction of J. G. Tappert. One of these was an original sugges- 
tion of my associate K. Levenberg.? The second method is the subject of this note. 

This method is not new. Levenberg found it set forth in a paper by Cauchy dated 
1847.* That it has become a standard procedure in analysis is clear from a recent paper 
by Courant.‘ Nevertheless it does not appear to be well known to authorities on nu- 


* Received Jan. 22, 1944. 

1 See, for example, W. E. Deming, Some notes on least squares, U. S. Dept. of Agriculture Graduate 
School, 1938, p. 31 ff., or E. T. Whittaker and G. Robinson, The calculus of observations, Blackie and 
Son, London, 1940, p. 214. Deming’s treatment is also given in his book, Statistical adjustment of data, 
John Wiley & Sons, New York, 1943, p. 52 ff. 

2 K. Levenberg, A method for the solution of certain non-linear problems in least squares, Quarterly of 
Applied Mathematics, 2, 164 (1944). ‘ 

3A. L. Cauchy, Méthode générale pour la résolution des systémes d’équations simultanées, Comptes 
rendus, Ac. Sci. Paris, 25, 536-538 (1847). 

‘R. Courant, Variational methods for the solution of problems of equilibrium and vibrations, Bull. 
Amer. Math. Soc. 49, 1-23 (1943). See especially pp. 17-20. Courant calls the method the “method of 
gradients” and ascribes its origin to a paper published by Hadamard in 1907, 





HASKELL B. CURRY 259 


merical computation. It was used for the case of linear equations by G. Temple;$ 
but he gave no reference to Cauchy’s work nor indeed to any previous use of the 
method. Accordingly there is room for an exposition of the method with emphasis on 
its practical aspects. 

This note also contains an outline of a convergence proof. Cauchy stated that the 
process converged but gave no proof, at least in the paper cited.. Temple’s conver- 
gence proof applied only to the linear case. Courant (l.c.) gives references to papers 
dealing with the method; but some of these were not accessible to me under wartime 
conditions. The convergence proof, as outlined here, is an elementary one and gives 
a weak result. 

The argument is, incidentally, capable of generalization to certain cases where 
there are infinitely many parameters, i.e., where G is a function of a vector x belonging 
to a suitable abstract space. The essential point is that there be a vector function 
H(x), the gradient, such that for vectors x, y and scalar ¢ 


d 
, [G(x + ty)] = H(x + ty)-y, 


where the dot indicates a scalar product.* Such generalizations will not be considered 
explicitly. 

2. Explanation of the method. The letters x, z will be used to stand for the n-tuples 
(vectors) (x1,- +--+, %n,) and (2, ---, Zn) respectively. It will be convenient also to 
think of the vector x as a point and z as a set of direction numbers of a direction, viz. 
the direction z, emanating from x. Superscripts will be used systematically to dis- 
tinguish different points and their corresponding directions. 

Let us suppose, then, that we start at a point x° and determine the direction in 
which G decreases most rapidly. This direction is given by 2; = —\ 0G/0x; or, in vec- 
tor form, z°= —X grad G, where } is an arbitrary positive factor of proportionality. 
(In practice we should either take \=1 or choose \ so that the vector z is of unit 
length.) Then the function g(#) =G(x°+éz°) has a negative derivative at ¢=0. It will 
therefore be possible to find a ¢>0 such that 


g(t) < g(0). (4) 


With such a ¢ we can take x! =x°+/z° as a new starting point and continue. We should 
then have a sequence of points x®, x}, x*, - - - such that G(x*t+!) <G(x*). Under suit- 
able restrictions (to be considered later) the sequence will attain or converge to a sta- 
tionary point of G. 

The determination of ¢ can be accomplished by trial. If no other indication is avail- 
ble we can take as first trial value the intercept of the tangent to the curve y =g(#) 
on the /-axis.’ If this fails to satisfy (4), it is too large; we can then take half of it, and 
so on. In this process we can draw a rough graph of g(t), and after a few trials it is 


5G. Temple, The general theory of relaxation methods applied to linear systems, Proc. Roy. Soc. 
London (A) 169, 476-500 (1939). For this reference I am indebted to H. Hotelling. 

* Even the case where such a gradient does not exist, there being only a total differential, can prob- 
ably be handled by a method which bears the same relation to the present method that Temple's method 
for gyrostatic systems does to his method of steepest descent. 

7 This was done by Cauchy (l.c.). It represents the approximation by Newton's method to the 
smallest positive zero of g(t). This is a reasonable guess for a G given by (2), where the numerical value is 
zero. In the least square cases (where G20) the guess is often many times too large. 











260 NOTES [Vol. II, No. 3 
usually possible to locate a ¢ which is at or near a minimum of g(t). Experience will 
presumably disclose many ways to shorten the process in individual cases. 

If we take for ¢ precisely the smallest positive root of 


g'(t) = 0, (5) 


the process has the following geometrical interpretation. Starting at x°, we determine 
the direction in which the surface 
y = G{x) = G(x, +++, Xn) (6) 


is descending most rapidly. We continue in that direction until we find ourselves 
going along a contour (i.e. a horizontal section of the surface). Then we stop, take a 
new direction of steepest descent, and so continue. Since the direction of steepest 
descent is always normal to the contour it follows that the directions z* and z*t! are 
at right angles.* This is important in the convergence proof. 

3. Proof of Convergence. Let us suppose now that G(x,-- - 
has continuous first partial derivatives at all points within or on the boundary of a 
region S. Let x° be a point within S. Let C be the broken line path starting at x° and 
going in the direction of steepest descent at x°® until it reaches either the boundary 
of S or the next approximation x! determined as in §3 with ¢ the least positive root 
of (5); in the latter case the broken line goes in the direction of steepest descent at x! 
until it reaches the boundary of S or the next approximation x? determined in the 
same way; and so on. Then G is monotone decreasing along C. There are three possi- 
bilities: (1) The path C may run into the boundary of S. (2) The path C may termi- 
nate at a point where the direction of steepest descent does not exist, i.e. at a sta- 
tionary point of G. (3) The process may continue indefinitely. The first possibility 
will certainly be excluded if the value of G at x® is less than at any point on the 
boundary of S. If the second possibility occurs the case is trivial. I shall make the 
limitation just stated in regard to G(x°) and shall suppose that the process continues 
indefinitely. 

Under these presuppositions let x* be a limit point of x®, x’, - 


that 


, X,) is defined and 


- +. Then it is clear 


G(x”) < G(x*) (k= 0,1,2,---). (7) 
It will now be shown that x” is a stationary point of G. 
Let us suppose the contrary. We write H(x) =grad G, h(x) = | H(x)| , and let z(x) 


be a unit vector; thus 
H(x) = — A(x)z(x). (8) 


According to the supposition, 4(x*) #0. Hence it will be possible to find a spherical 
neighborhood U of x* such that for x in U, 

| H(x) — H(x*) | < eh(x*). 
Then it will follow that | h(x) —h(x*) | <eh(x*), |z—z"| <2e. Hence, from the fact 
that z* and z*+! are at right angles, one can conclude that if x* isin U, x**! is certainly 


not in U (provided « is not too large). 
Next, let K be the conical sector of U for which 





8 This is also easily proved analytically. 
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Then it may be shown, by reasonably straight-forward methods, that for x in K, 
G(x) < G(x”). (9) 


Let V be a subneighborhood of U such that for x in V the ray in the direction z (as 
given by (8)) from x intersects K. Such a V exists if ¢<4. Since x” is a limit point, 
there exists an x" in V. Then the ray in the direction z" from x* will have a point y 
in K. This y cannot be beyond x**!, since x"*! is not in U and U is convex. Hence, 
by the monotonic character of G on C and by (9), G(x*t!) <G(y) <G(x*), which con- 
tradicts (7). This contradiction came from the assumption that h(x”) #0. 

The following example shows that we cannot expect a better result without further 
restrictions on G. Let G(x, y) =0 on the unit circle and G(x, y) >0 elsewhere. Outside 
the unit circle let the surface have a spiral gully making infinitely many turns about 
the circle. Then the path C will evidently follow the gully and have all points of the 
unit circle as limit points. 

In a practical problem however we often know in advance that there is a unique 
minimum of G within S; in these cases convergence is assured. If G is given by (2) 
and the Jacobian of the f’s does not vanish in S, then every stationary point of G 
is a solution of (1); if there is only one such solution the process converges to it. 

4. Concluding Remarks. In regard to the practical aspects of the method the fol- 
lowing points are to be noted: (1) It does not require any calculation of second deriva- 
tives. This is important for the application mentioned, where these second derivatives 
are numerous and complicated. (2) It involves only direct calculations of G and its 
first derivatives. (3) The approach to the limit, if any, is along a path C consisting of 
straight line segments, adjacent segments being approximately at right angles. 

A comparison with Levenberg’s method in regard to these three respects is now in 
order. In the first respect the two methods are alike. In the second respect Levenberg’s 
method is more complicated because each stage requires the solution of a set of “nor- 
mal” equations, as in the traditional method of least squares. In the third respect 
Levenberg’s method is like the present one in that it involves approach along a broken 
path C; but the individual pieces of C are curved. There is evidence that in practical 
problems these curves follow the natural valleys of the surface, so that each step 
brings us further toward the goal. As to which of these opposing characteristics is 
the more important is not yet settled.® 

Another point is that the process is not invariant under certain elementary trans- 
formations, e.g., a change of scale of one or more of the x;. Thus, if the surface (6) (for 
n =2) is a hemispherical bowl, the direction of steepest descent is along the meridian, 
and the minimum is reached in one step. If, however, the scale on the x;-axis is 
changed so that the surface becomes ellipsoidal, the direction of steepest descent is 
no longer directed toward the minimum. The suggestions inherent in this lack of in- 
variance have not yet been fully worked out. 


® The engineers at Frankford Arsenal prefer the Levenberg method for the problems which have 
confronted them; but I do not know to what extent they have exploited the present method. Since the 
Levenberg method is confined to a G of the form (3) and makes use of that representation, it would not 
be surprising if it should prove superior for that case. On the other hand it is not difficult to concoct arti- 
ficial examples for which the method of steepest descent is superior in the third respect as well as the sec- 
ond, at least for certain determinations of the weighting factors. 
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ON THE BENDING OF A CLAMPED PLATE* 
By A. WEINSTEIN (University of Toronto) and D. H. ROCK (Rhode Island State College) 


The present paper contains an application of a recently developed variational 
method! to the boundary value problem of the bending of a clamped plate of arbi- 
trary shape. It will be shown that this problem can be linked to the simpler problem 
of the equilibrium of a membrane by a chain of intermediate problems, which can be 
solved explicitly and in finite form in terms of the membrane problem. In the inter- 
mediate problems, the deflection converges uniformly in the domain of the plate 
(including the boundary) to the deflection of the clamped plate, and the derivatives 
of all orders of the deflection converge uniformly in every domain completely interior 
to the plate. (In the Ritz method, not even the convergence of the slopes can be 
guaranteed.”) The method yields numerical results for plates of all shapes for which 
the membrane problem (which we shall call the base problem) admits an explicit solu- 
tion. As an example we shall consider a clamped square plate under a uniform load. 
This problem has been the object of numerous investigations,’ some of which are 
theoretical, while others are purely numerical, use infinite simple and double series, 
and operate with an infinite number of linear equations and an infinite number of 
unknowns.‘ An inspection of the general formulae derived in the present paper, for- 
mulae which become simple in numerical applications, would show how some of the 
numerical methods might be rendered rigorous. The convergence of higher deriva- 
tives is of great practical interest for the approximate computation of the stresses. 
(Cf. Handbuch der Physik, Springer, Berlin, Vol. VI, 1928, pp. 220-221.) 

Let us denote the domain of a clamped plate by S and its boundary by C. The de- 
flection w(x, y) corresponding to a load g(x, y) and to a flexural rigidity D satisfies 
the differential equation 

AAw = q/D (1) 


with the boundary conditions 


w= 0, (2) dw/dn = 0, (3) 


on C. It is well known that w is the solution of the variational problem VP, 


J(w) = ff [(Aw)? — 2gw]dxdy = min., 
s 


* Received Feb. 8, 1944. 

1 A. Weinstein, Mémorial des Sciences Mathématiques, No. 88, 1937. 

2 Cf. R. Courant, Variational methods for the solutions of problems of equilibrium and vibrations, Bull. 
Amer. Math. Soc., 49, 1-23, 1943, especially p. 11. See also K. Friedrichs, Math. Annalen, 98, 217, 1928. 
The method of finite differences does not give satisfactory numerical results for clamped plates. 

3 The bibliography given in S. Timoshenko, Plates and shells, McGraw-Hill, 1940, p. 222 covers 
papers from 1902 to 1939 and shows the persistent interest in this problem. 

4H. W. March, Trans. Amer. Math. Soc., 27, 307-317, 1925, proves the weak convergence of the 
approximations given by his series. Cf. I. S. Sokolnikoff, Mathematical theory of elasticity, Brown Univer- 
sity, 1941, p. 387. 

5 Cf. for instance a note by C. Miranda, Rend. Semin. Mat. di Roma, 1, 262-266, 1937. 
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where ¢=q/D and the boundary conditions (2) and (3) are in effect. By withdrawing 
the condition dw/dn =0, we obtain the variational base problem VP 9: J(wo) =min., 
with the condition wo=0 on C. The corresponding differential equations problem 
DP, is: AAwy =, with the boundary conditions w)=0 and Aw,)=0 on C, the latter 
condition being a natural boundary condition which is automatically satisfied by a 
solution of VP». It is well known that DP, can be solved in terms of the problem of 
the equilibrium of a membrane. In fact, putting Awo=fo in DP» (wo=0 on C), we 
have Afy=g@ in the domain S and fy) =0 on C.§ Let us denote the solution of the equa- 
tion Aw ) =f) with the boundary condition wo=0 by wo=Gfo. Thus we have fo =Gq@, 
whence wo = GG@. (The formula wy = Gf, can be written explicitly in the form 


are f J _8( 4 WB Wfolbs sdb, 


where g(x, y, £, 7) is the Green’s function for the domain S.) 

We now link VP, to VP by a chain of intermediate variational problems 
VP, VP2, +--+ introduced in the following way. We let ;(x, y), po(x, y),-°* de- 
note a complete (bit not necessarily orthogonal) sequence of linearly independent 
harmonic functions in S. (It has been shown, I.c.,!1 how a sequence of this kind can 
be derived from the solutions of the problem of a vibrating membrane for any do- 
main S.) 

VP,, (m=1, 2,-+- +) is then defined as the problem of finding the solution w,, of 
J(wm) =min., with the boundary conditions 





: dWm 

Wm =O0onC, Pr ds = 0, k=i,2,---,™. 
c dn 

By Green’s formula these conditions can be replaced by the conditions 


Wm = 0 on C, (pr, AWm) = 0, k=i,2,---,™, (4) 


(Pry Awm) = ff prAwndxdy. 
Ss 


(Similar notations for the “scalar product” of two function, like p, and Awm, will be 
used throughout this paper.) The corresponding differential problem DP,, is given by 
the differential equation 


where 


AAwm = @ (5) 

with the conditions (4) and the natural boundary condition 
AWm = Gmipi + Gmep2 + +++ + dmmpm on C, (6) 
where Gm1, @m2, °° * » mm are constants to be determined. The solutions of DP,, can 


be easily obtained in terms of solutions of the membrane problem already used to 
solve DP». In fact, putting Aw», =fn (Wm =0 on C), we have, in our notation, Ww», =Gfn. 


Also, we obtain from (5), 
Afm = q; (7) 





6 For rectangular plates DPo is the problem of the supported plate. 
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and from the boundary condition (6), 


ie = Qmip1 + Im2p2 + aa + AmmPp m- 


(8) 


These can be written as follows: 


a(t -> an:ps) =ginS (7’) 


im 
In — > Amipi = 0 on C. (8’) 

Therefore we have in S sas 
fm — > Imipi = Gq, (9) 

inl 


and since Wm = Gfm, it follows that 


Wm = GGG + >> amiGPi, (10) 


i=! 


where GGG =» is the solution of the base problem. 
The conditions (pz, Aw») =0 yield, in view of (9), the following system of m linear 


equations for the m constants @m1, @m2, ** * , @mm: 
m 
Lo ami(pis Px) = — (G,Gpr), k= 1,2,-++,m, (11) 
i=1 


which can be solved, since their determinant is Gram’s determinant of the independ- 
ent functions pi, p2, - + + , Pm, and is different from zero. 

In another paper, based on a previously developed method for the computation 
of frequencies and buckling loads,’ it will be proved that the approximate solutions w», 
and their first derivatives converge to the deflection and slopes of the clamped plate. 
Here we shall apply our formulae to the case of a uniformly loaded square plate. 
The domain S will be defined by the inequalities |x| <7/2, | y| S7/2. We put 
g=q/D=1. Since the deflection of a uniformly loaded square plate is symmetrical 
with respect to the coordinate axes, we may use a sequence of even harmonic func- 
tions p;(x, y) as given by (12) below.* All computations can be performed without 
the use of Green’s function for the square. The deflection wo of the supported plate 
is given by the well known formulae of Navier. 

Calculation of wm for the uniformly loaded square plate. 


We use the set of functions 


cosh a;x cos a;y + cos a;x cosh ayy : 
Sh aX COS a;4 S a ay apes ee (12) 





pi(x, y) = 
cosh (a;7/2) 


7N. Aronszajn and A. Weinstein, On the unified theory of eigen-values of plates and membranes, 


Amer. J. Math., 64, 625-645, 1942. 


8 For non-uniform loading, the sequence given l.c. 1, {16 must be used. 
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and denote Gp; by 2;. Then, by the definition of G, we have Av; =AGp; = p; in S, with 
v;=0o0n C. 
If we set 0; =u,+4%2, where 


cos ay COS ajx 


*, = X(x) - ——_—_» 2 = ¥(y) ————_. 
cosh (a,7/2) cosh (a,7/2) 





then X(x) and Y(y) satisfy the differential equations 
X” — aX = cosh a,x, Y” — a?Y = cosh ay, 
with the boundary conditions 
X(+ 2/2) = 0, Y(+ 2/2) = 0. 
The general solution for X (x) is 
X(x) = 


ra x sinh a;x + A cosh a;x + B sinh a;x, 
4Q; 


where A and B are determined by the boundary conditions. Hence 


1 


Tv 
xX =— | sinh ax — ry tanh }a;7 cosh a | 


Qa; 
The solution for Y(y) is given by a similar expression, so finally 


. J 


Tv 
—-- ——-- E sinh a;x — — tanh $a, cosh a | COS ay 
2a; cosh (a, /2) \ 2 


T 
+ E sinh ayy — > tanh 3a,;r cosh ay | cos aa} . (13) 


Using (13), we obtain the general formula for (g, Gpi) = (1, v3), 





4 sin (a;r/2) Px 1 
(1, »;) = —_——| — sech? }a;r — — tanh ja,r |. (14) 
a’ 2 Qi 
For (pi, px) we have 
Sarsouz(—) *+* 
anne, teh 
(a; + ax)? 
T 1 2 F 
= |= sech? jayr + — tanh fair | + —> ¢= k. (15) 
2 Qi a; 
From (14) we find that 
(1, 21) = — 2.670644, (1, v2) = 0.049299, (1, v3) = — 0.006400, 


(1, v4) = 0.001666, etc., 


and from (15) we get the following table for (p;, pe), (¢, R=1, 2, 3, 4): 
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NA 1 2 3 4 
1 5.665118 | —0.240000 0.059172 | —0.022400 
2 1.270836 | —0.103806 | 0.049941 
3 | 0.708321 | —0.051132 
4 | | 0.489615 


The equations for the determination of a,,; are then 
5.665118a;, = 2.670644, for ay, 
5.665118¢2; — 0.240000a2. = 2.670644 
for do, and a 
— 0.240000a2; + 1.270836a22 = — 0.049299 


22) 


These yield the successive values 


a31 > 0.471419, 


ao, = 0.473564, doo = 0.050641, 
a3; = 0.473728, a3. = 0.048762, a33 = — 0.023392, 
a4, = 0.473749, dg. = 0.048394, a43 = — 0.022656, a4, = 0.010970, 


Then, since the solution for the simply supported plate is given by 


: 16 sin m(x + 2/2) sin n(y + 2/2) 
w(x, vy) =—)>, >, \ /2) hein ’ (m,n = 1,3,5,---), 
oe mn(m? + n*)?- 





the successive approximations to the deflection are: 
WwW; = Wo +- 0.4714192,, 
wo + 0.4735640; + 0.05064100, 


We 


The maximum deflection, which occurs at the center of the plate, is found to be 
0.123342 when we is used. The next approximation affects only the fourth significant 
figure. 

A calculation of the normal derivative of w», (m=0, 1, 2, 3, 4) along x =72/2 for 
values of y from 0 to 7/2 at intervals of 1/16 yields: 





m| y=0 r/16 24/16 3x/16 4n/16 54/16 | 6/16 | 7/16 | /2 
0 | —.41795 | —.41087 | —.38954 | —.35409 | —.30519 | —.24382 | —.17095 | —.08839 | 0 
1 | —.00763 | —.00686 | —.00445 | —.00062 | +.00443 | +.00912 | +.01221 | +.01100 | 0 
2 | +.00197 | +.00122 | —.00047 | —.00189 | —.00208 | — .00060 | +.00216 | +.00355 | 0 
3} —.00053 | —.00019 | +.00044 | +.00055 | —.00020 | —.00093 | —.00018 | +.00151 | 0 
4 | +.00029 | —.00003 | —.00031 | +.00008 | +.00036 | —.00021 | —.00033 | +.00066 | 0 


The maximum value of the slope in the interior of the plate is found to be about 


—0.122; this occurs at x =52/16, y=0. A comparison with the maximum deviation 
of the normal derivative along the edge for m =4 shows that the latter is less than 1% 


of the maximum slope in the interior of the plate. 
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THE DYNAMICS OF A DIFFUSING GAS* 
By HENRI PUTMAN (Université Laval, Québec) 


By use of a hydrodynamical approach, Stefan! derived the following equation 
for the diffusion of two gases: 


Op. 
piti = ———A 12 1P2( M1 _ Ue). (1) 
Ox 


Here /; is the partial pressure of the first gas, p; is its density, m, £; are respectively 
the x-components of the velocity and acceleration of one of its particles, uz, pz refer 
to the second gas, and Aj is a constant. There are two other equations similar to (1) 
corresponding to the y and z-components, and a further set of three equations for 
the second gas. 

Equation (1) is a simplified form of Maxwell’s equation of diffusion.? It states that 
there acts on a particle of the first gas a force due to the pressure gradient of the first 
gas, and a force proportional to the difference of the velocities of the two gases. 

The ordinary equation of diffusion, or Fick’s law,* was deduced by Stefan! from 
(1) and the equation of continuity by assuming that & was negligible. 

We shall now assume that m% is negligible. This is the case in which the second 
gas is immobile and the first gas diffuses through it. Some problems involving two 
gases can be reduced to just such a problem.‘ 

If now in (1) we set u,=v, £;=dv/dt, Aype=a, pi=p, pi=p, we obtain 


dv 1 dp 


The corresponding three dimensional form when there is present in addition a body 
force per unit mass represented by the vector F is 


d 
edie i ee ea (2) 


dt 7 p 
where v is the velocity vector. This equation is the equation of motion of a viscous 
fluid with the viscosity terms replaced by a force proportional to the velocity. 
We shall now deduce some additional equations which are consequences of (2). 
If we multiply (2) scalarly be an arbitrary virtual displacement de, we obtain, 
since 6=V-ée, 
dv bp 
—-d5e = F-5de — — — av-ée. (3) 
dt p 


* Received Feb. 8, 1944. 

1 J. Stefan, Ber. der Wiener Akad. 63 (2), 63-124 (1871). 

2 J. C. Maxwell, Phil. Mag. (4), 35, 185-217 (1868). 

*R. M. Barrer, Diffusion in and through solids, Cambridge University Press, Cambridge, 1941, p. 1. 

‘ B. Lewis and G. v. Elbe, Combustion flames anc explosion of gases, Cambridge University Press, 
Cambridge, 1938, p. 224. 
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When the expression v-d(de) /dt = 6(4v?) is added to both sides of (3), we obtain 


d bp 
= (v-de) = F-de — — + 5(40?) — av-ée. 
at p 


(4) 


If now we set de =de, where de follows the natural motion of the system, and assume 
that F has a potential U, then from (3) 





dp 
d(3v?) = dU — — — av-de, (5) 
p 
or 
dp 
d($v”) + av*dt = dU — — - (6) 
p 
If the temperature of the gas is constant, we have 
dp 
p = Kp, ft-xime, 
p 
K being a constant. Integration of (6) then yields 
= const. (7) 


3° + Kinpt+ af v-de 


When the only body force is due to gravity, and the x-axis is vertically downward, 


U=gx and (7) becomes 
g iy 


K v* a 
= ~ f v-de = 
g a 


i~—he~— = 


const. 


(8) 


We now return to (3). If there is a straight or curved axis of symmetry such that, 
if s is the arc length of this axis, v, p, p are functions of s, ¢ only, and if we set de =u,és, 
where u, is a unit vector tangent to the axis of symmetry, then (3) becomes 


dv op 
—-u,ds = gix — — — av-ut,és, 
dt p 
or, since Vv: U, =2, 
dv Ox 1 dp 
—=s¢—=—-— — — av. (9) 
dt Os p os 
Now 
Ov Ov dv Ov Ov 
dv = — ds + — di, —= — y+ — 
Os ot at Os ot 
whence (9) takes the form 
Ox 1 dp 90 fv 1 dv a 
See th Sat lk lysate ares ls sank es le tii (10) 
Os pg Os Os \2g g ot g 
The equation of continuity is 
Op 0 
2 — + — (pQv) = 0, (11) 


ot OS 
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where @ is the area of the cross section. When the gas is at constant temperature, 
pb =Kop, and (10) and (11) become 


°(2-=inp)==(=)4+— 54", (12) 
Os g Os \2g g ot g 
ae ba (pQr) = 0. (13) 
Ot Os 


If the flow is steady, 0v/dt=0 and (12) becomes 


a v v a 4 
n-mt—in = S44 f ods, (14) 
g§ pi 2g gw gd 
where the subscripts zero and one refer to two cross sections, both viewed at the 
same instant. We set P=3(po+)p1), T=po— pi, a=}3m/P, whence 


po = P+ })x = P(1 +a), pi = P — }x = P(1 — a), 
Po li+a 
In— = 


In 

pi 1 — 
When the difference of pressure at these two cross sections is small, $a? is much 
smaller than 1 and we have approximately In (po/p1) = 2a=2/P. If w is the specific 


weight at pressure P, then 





= 2(a+ fa*®+---) = 2a(1+ ja? +---). 





Kw K 1 K _ . 
P = Kp = —> > i a= * Si ee 
g gp 2 26w g pi w 
and (14) becomes 
2 2 
_ v v a a 
as age re | vds. (15) 
w 2g 2g gw 


Let us return once more to (2), and operate on it with VX and V.-, to obtain 


d fa @ a 
“(*)=(*-4)r-<« (16) 
dt\p p p 


dé 
= = Ay — a0 — 1° [(vi-V)v], (17) 
¢ 


where o=V Xv, 6=V-v, Y= U—p/p, p is assumed to be a function of p only, and 
ov OV 
Vi: [(vi-V)v] = @ — 230k: (= xX =), 
Ox dy 
i, j, k being unit vectors along the x, y, z-axes, respectively. The equation of con- 


tinuity is 


dp 
dt : 
or 
Op 
— = — V-(pv). (19) 


ot 
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If a is not a constant, but a function of e(x, y, 2) for example, we must add the 
term —(vXVa)/p to the second term in (16). 

We shall now examine the propagation of discontinuities in the boundary condi- 
tions, in the simple case of one dimension without gravity. Equations (2) and (18) 
become 


ies tuned (20) 
—_—- — a av = ’ 
ot Ox p Ox 
Op Op ov 
—+—v+p—=0. (21) 
Ot Ox Ox 


Ov Ov K op 
—+y7—+— — = - 4, (22) 
ot Ox p Ox 
Op Op Ov 
—+7—+ p— =), (23) 
ot Ox Ox 


which are a system of two simultaneous quasilinear equations of the first order. These 
equations are usually referred to as Hamburger’s equations in two dependent varia- 
bles. By standard procedures,‘ the ordinary differential equations for the characteris- 
tics are found in the form 


dt 


- Fa = (vo+ K!/2)-!, dv + K'!/*2p-'dp = — av(v + K"/?)—dx, (24) 
dx 
dt 
rs = (v — K?/?)-1, dv — K*/2p9—'dp = — av(v — K"/?)—dx. (25) 
ax 


If x, ¢, v are regarded as a rectangular cartesian coordinate system, and x, ¢, p as 
a second system, the solutions of Eqs. (24) and (25) can be represented graphically 
as surfaces. We may assign boundary conditions as follows: 

(a) In both the xv- and the xp-planes, a curve is given. 

(b) In the ép-plane, a curve is given; the corresponding curve in the /v-plane must 
be determined by means of the characteristics. 

(c) In a plane x =const. =d of the xtp-system, a curve is given; the corresponding 
curve in the plane x =d of the xvt-system must be determined by means of the charac- 
teristics. 

Let us consider, for example, a gas for which K = 13.7 X 10‘m.’*sec.~?, a = 2.75 K 10° 
sec.—!, with the following boundary conditions: 

(a) In the xv- and xp-plane, v =v) =0, p = po = 0.073 kg.m.—‘sec.? This corresponds 
to an initial pressure at rest of 10000 kg.m.~? 

(b) In the ¢p-plane, the curve p =f(t) is given. 

(c) In the plane x = 2000 m. of the xtp-system, the curve p=0 is given. 

When v=0, we find from the first of Eqs. (24) that dt/dx =K-/?=370 m.sec.—'. 
Thus the time that the initial discontinuity requires to cover the distance x = 2000 m. 
is 5.4 sec. 


5 A. R. Forsyth, Theory of differential equations, Cambridge University Press, 1906, vol. 5, p. 435. 
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ON THE TREATMENT OF DISCONTINUITIES IN 
BEAM DEFLECTION PROBLEMS* 


By ERYK KOSKO (Ecole Polytechnique, Montréal) 


In a note on the treatment of discontinuities in beam deflection problems (Quar- 
terly of Applied Mathematics, 1, 349-351) Mr. C. L. Brown suggests the use of 
Heaviside’s unit step function. He thus avoids what he calls the “sectionalizing” 
treatment in the integration of the differential equation for the deflection of a beam 
with discontinuous transversal loading. 

Mr. Brown's method appears to be equivalent to the procedure which seems to 
have been first developed by R. Macaulay! and has since been included in several 
British textbooks.?* In order to establish expressions for moments with discontin- 
uous variations, Macaulay introduces terms in twisted brackets, such as {x—a}, with 
the convention that these terms be neglected when the quantity within the brackets 
becomes negative. When integrating the term in question, the quantity in brackets 
is to be regarded as the independent variable instead of x; the indefinite integral of 
{x—a} would be }{x—a}?. 

Taking Mr. Brown's example, the expression for the bending moment of the beam 
(Il. c., Fig. 1, p. 349) would be with Macaulay’s notation: 

d*y 
ai = M = —Mi+ Rx — P{x—- a}. 
x? 





The first integration would give 


ly 
El — = — Myx + 4Rit* — 4P{ x — 0}? +C, 


dx 





and the second integration 
Ely = — 3M,x* ARix* — 2P{x — a}? +Cyx+ Cr 


All the above equations hold at all parts of the span so that there are only two 
constants of integration to be determined from the conditions at both ends of the 
span, instead of having two for each section of the span as in the classical treatment. 

It is therefore apparent that Macaulay’s twisted bracket is but another symbol 
for the multiplication by the unit step function. 

An important remark is to be made about the use of this procedure with regard 
to distributed loads. These must always be made to extend to the right-hand extrem- 
ity of the beam, introducing negative loads if necessary. An extension of the method 


* Received March 27, 1944. 

1 R. Macauly, Note on the deflection of beams, Messenger of Mathematics, 48, 129 (1919). 
2 R. V. Southwell, An introduction to the theory of elasticity, Oxford, 1936, §§194-196. 

3 J. Case, The strength of materials, 2nd edition, Arnold & Co., London, 1932, §169. 
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due to H. A. Webb’ covers the effect of a concentrated bending couple applied at an 
intermediate point of the span. As an application, R. V. Southwell suggests (Ex- 
ample 14, I.c.) the derivation of the theorem of three moments for a continuous beam 
by the same method. The use of the method for a beam with a stepwise variation of 
bending rigidity seems however to be Mr. Brown’s original contribution. 

It is hoped that the discussed method, whatever the symbols used, will get more 
attention from engineers on this side of the Atlantic, as it carries with it a very sub- 
stantial shortening of the computations. 


FORMULAS FOR COMPLEX INTERPOLATION* 
By A. N. LOWAN anp H. E. SALZER, (Math. Tables Project, Nat. Bureau of Standards) 


An analytic function of z=x+7y may be approximated by a complex polynomial 
of degree m passing through +1 points in accordance with the Lagrange-Hermite 
formula of interpolation. For the important special case when the given +1 points 
are equidistantly spaced along any straight line in the z-plane, the following tables 
give the real and imaginary parts of the coefficients A;(P) of the interpolation 
polynomial f(z) =Axz(P)f(z:), where P = (z—20)/h=p+ig and h is the complex tabular 
interval. The formulas cover the cases ranging from complex quadratic (3 points) 
to complex quintic interpolation (6 points). 


Quadratic interpolation (3 points) 


ReA_,(P) = 3[p(p — 1) — q’], 3mA_,(P) = q(p — .5), 
ReA(P) = 1— p?+ 9g’, ImAo(P) = — 2pq, 
ReA(P) = 3[p(p + 1) — q?], 3mA,(P) = q(p + .5). 


Cubic interpolation (4 points) 


(1-9) 
ReA_(P) = —f [p(p — 2) — 3q°], 3mA_,(P) = < [g? —2+3p(2 — p)], 
) ) 
ReAo(P) = 1+4[p(p? — 2p — 1) + g°(2 — 3p)], ImAd(P) = = [p(3p — 4) — q?—1], 
ReAi(P) = — 3[p(p — 2)(p +1) + 9°(1 — 3p)], ImAi(P) = J [p(2 — 3p) +94?+2], 
y eH (gs — 3g a 
ReA(P) = 4 (p? — 3g? — 1), 5mAx(P) = = fiat gt — §] 
) 


* Received April 17, 1944, 
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Quartic interpolation (5 points) 


1 
ReA_o(P) = [p(p? — 1)(p — 2) + g2(q? + 1) + 6pqX(1 — p)], 
3mA_o(P) = ze [q2(1 — 2p) + 2p? — 3p? — p+ 1], 


12 
ReA_(P) = — ; [p(p — 1)(p? — 4) + g°(q? + 4) + 3pq?(1 — 29)], 
jmA_,(P) = — < [4p3 — 3p? — 86 + 4 — q?(4p — 1)], 

ReAo(P) = ; [(p? — 1)(p? — 4) + g°(q? — 6p? + 5)], 

5mA,(P) = ” (2p? — 2g? — 5), 

ReA(P) = — - [p(p + 1)(p? — 4) + 9g? + 4) — 3pq°(1 + 29)], 
SmA\(P) = — — [4p* + 3p* — 8p — 4 —gX(1 + 49)], 

ReA(P) = 7 [p(p + 2)(p? — 1) + g°(q? + 1) — 6pg*(1 + 9)], 


3mA(P) = = [2p ++ 3p? — p—1 — (2p + 1)]. 


Quintic interpolation (6 points) 
1 
ReA_»(P) = co — p(p? — 1)(p — 3)(p — 2) + 5q%(p — 1)(2p? — 4p — 1) + SqX(1 — p)], 
SmA_a(P) = —— [— 5p(p — 2)(6* — 2p — 1) — (g" — 6)(Q" + 1) + 10pa%( — 2)], 
, p : , q 
ReA_(P) = ae — 1)(p? — 4)(p — 3) — g*(10p? — 24p — Sq’ — 3)| - ras + 4), 
SmA_4(P) = - [(g? + 4)(q? — 3) + p(Sp* — 16p? — 3p + 32) + pg?(16 — 109)], 


1 
ReAo(P) = — [(p* — 1)(p* — 4)(3 — p) + g*(h0p* — 189" — 15p + 15) + a3 — 5p)], 


SmAo(P) = — = [(g? + 1)(g2 + 4) + p(Sp* — 12p? — 15p + 30) + pg*(12 — 109)], 
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ReA(P) = - [p(p + 1)(p? — 4)(p — 3) + g2(— 10p* + 129? + 21p — 8) + g4(5p — 2)], 
5mA,(P) = rn [(o2 + 4)(¢2 + 3) + p(Sp* — 8p? — 21p + 16) + pg%(8 — 109)], 
ReAx(P) = 7 [— p(p?— 1)(p— 3)(p + 2) — g°(— 10p° + 6p? + 21p — 1) + (1 — 5p) ], 
3mA,(P) = — mr [Spt — 4p8 — 21p? + 2p + 6 + g2(g? + 7) + 2pg°(2 — 5p)], 


ReA3(P) = oe [(p? — 4)(p? — 1) + 5q°(q? — 2p? + 3)], 


Il 


— [(g? + 4)(q? + 1) + 5p°(p? — 2g? — 3)]. 


120 
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BOOK REVIEWS 


Ten lectures on theoretical rheology. By Markus Reiner. Nordemann Publishing Co., 
Inc. New York, 1943. iv+164 pp. $4.50. 


It is customary to define Rheology as the science of flow and deformation of matter. If this definition 
is taken literally, Rheology becomes practically identical with Mechanics of Continua, and it becomes 
hard to understand why a new term has been coined. It seems that a more adequate definition would state 
Rheology to be a Mechanics of Continua in which the ideal elastic body and the perfect fluid are almost 
as systematically disregarded, as they are over-emphasized in classical Mechanics of Continua. This 
undue prominence which the classical scheme gives to the ideal elastic body and the perfect fluid tends 
to be reflected in textbooks of Mechanics as well as in engineering curricula. The necessity of developing 
and studying the mechanics of other solids and fluids must therefore be stressed, even to the extent of 
creating a new term for this part of Mechanics of Continua which, up to a fairly recent past, has been so 
badly neglected. 

There is a definite need for a treatise covering the entire field of Rheology rather than parts, such as 
plasticity or dynamics of (Newtonian) viscous fluids. The present book is in the nature of an introduction 
to such a treatise. Four chapters (1-3, 9) deal with the analysis of stress and strain and the important 
decomposition of the tensors of stress and strain into isotropic and deviatoric parts. The author bows to 
convention in defining the strains as ¢z2=0u/Ox, +++ , €zy=0v/dx+0u/dy, +++. This procedure, a rem- 
nant from the time when the tensor character of strain and its geometrical implications was not yet fully 
realized, makes it necessary to denote the components of the strain tensor by ézz,+++ , }é@zy,*** , and 
deprives many relations of their natural symmetry. Four further chapters (4, 6, 8, 10) are devoted to. the 
discussion of certain rheological idealizations, viz. the ideal elastic solid, the (Newtonian) viscous fluid, 
and the materials customarily named after Maxwell (viscous fluid with relaxation of stresses), Voigt 
(visco-elastic solid), Saint-Venant (perfectly plastic solid), and Bingham (viscous fluid with yield limit). 
The remaining two chapters (5, 7) are concerned with the solution of special problems (tension and simple 
flexure of a prismatical bar), Einstein’s law of the viscosity of sols, and rheological models. 

Theoretical Rheology is a subject which cannot be treated satisfactorily without using the tool of 
tensor analysis. The author did obviously not want to suppose the reader to be familiar with this tool. 
On the other hand, limitations originally imposed on the time available for his lectures seem to have pre- 
vented him from presenting even the basic conception of tensors ina precise form. One wonders whether, 
under these circumstances, it would not have been preferable to restrict the discussion to the mechanical 
behavior of solids and liquids in pure shear. As it is, the uninitiated reader cannot fail to get the impres- 
sion that any nine quantities neatly arranged between double vertical bars constitute a tensor. The laws 
of transformation are touched upon in chapter 9 only, and are nowhere stated to form an integral part 
of the definition of tensor. On the other hand, the reader familiar with tensor analysis might wish to have 
a tensorial expression of the yield condition of plastic materials which is more adequate than the cryptic 
relation po =0 (p. 111, Eq. (4)), where po denotes the stress deviator and # is defined as “the yield stress.” 
Many other instances could be cited where the clarity of exposition has obviously suffered from the tend- 
ency to cram too much material into a text which, according to the preface, is intended as a brief intro- 


duction. In spite of such occasional shortcomings the book, which fills a patent need, will prove very useful. 
W. PRAGER 


Table of reciprocals of the integers from 100,000 through 200,000. Prepared by the 
Mathematical Tables Project, Work Projects Administration of the Federal 
Works Agency; conducted under the sponsorship of the National Bureau of Stand- 
ards. Official Sponsor: Lyman J. Briggs; Technical Director: Arnold N. Lowan. 
Columbia University Press. New York, 1943, viii+201 pp. $4.00. 

These tables are a useful supplement to the existing tables of reciprocals. The tabular interval is 
small enough to permit linear interpolation throughout; the differences decrease slowly from 100 to 


25 units of the last place. The arrangement of the tables is very practical. Seven significant figures are 
given (if a number has & figures before the decimal point, its reciprocal has k—1 zeros after the decimal 
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point, before the first significant figure). Moreover, the tables indicate the direction in which the last digit 


is rounded; this practical device reduces the relative tabular error to 2.5 X10-8, 
W. FELLER 


Table of the Bessel functions Jo(z) and J,(z) for complex arguments. Prepared by 
the Mathematical Tables Project, Work Projects Administration of the Federal 
Works Agency; conducted under the sponsorship of the National Bureau of 
Standards. Official Sponsor: Lyman J. Briggs; Technical Director: Arnold N. 
Lowan. Columbia University Press. New York. 1943, xiv+403 pp. $5.00. 

Many problems of mathematical physics and mechanics lead to Bessel functions of various types. 
The most important of these problems are sketched in the foreword to the present tables, written by Pro- 
fessor H. Bateman of California Institute of Technology. The number of existing tables of Bessel func- 
tions is also legion. The present tables contain a valuable bibliography listing some 65 tables of Bessel 
functions of orders zero and one. However, the new tables are unique both in range and extent. 

The Bessel functions J,(z) are defined by 

o (- 1)*z”+2* 
Ie) =U 
kmo RIT (vy + Rk + 1)2”*2% 





_and satisfy the differential equation 
2*u''(z) + 2u’(s) + (2? — v*)u(z) = 0. 
One has the recurrence formula 
Jy41(2) nl y-1(2) + 2vJ,(z)/z 

which enables one to compute J,(z) for all integers v given the values of Jo(z) and J;(s). These are now 
tabulated to ten decimal places. The entries are written in the polar form z= p(cos ¢+ sin ¢). The func- 
tions are tabulated along the rays ¢=0°, 5°, 10°, - - - , 90° for values of p from 0 to 10 in the steps of .01. 
The values of the functions in all other quadrants follow easily by means of simple symmetry relations. 

To facilitate interpolation, the book contains also tables of the coefficients in Lagrange’s interpola- 
tion formula which uses five equally spaced points. The coefficients are tabulated to 10 decimal places, 
the argument varying in steps of .001. These tables will, of course, be useful for many computations and 


are quite independent of the main tables. 
W. FELLER 


The methodology of Pierre Duhem. Armand Lowinger. Columbia Univ. Press, 184 pp., 


1941. $2.25. 

The French theoretical physicist Duhem (1861-1916) devoted his life-work to thermodynamics, 
although he is mostly remembered today for his studies in medieval science. He moreover published his 
views on the methods, aims and significance of physics in a few scattered papers and in one connected 
account, a book entitled La Théorie physique, son objet et sa structure (1906). Mr. Lowinger has given an 
excellent presentation of Duhem’s ideas. These are challenging, for on the one hand Duhem was an 
eminently “classical” physicist, strongly opposed to atomism, opposed also to Maxwell’s electromagnet- 
ism, so that from our present point of vantage we can prove him wrong on both these counts; but on the 
other hand, he believed with Kirchhoff that physical theory is a description, not an explanation; with 
Mach, that its purpose is intellectual economy, and so was on the way, with these and other correlated 
ideas, towards present-day scientific pragmatism. Duhem’s opinions on physical method were rooted in 
his metaphysical beliefs, a fact obvious to his readers and critics, but which he was very anxious to deny. 
Mr. Lowinger’s presentation of Duhem is faithful and unbiased; he gives us his own views in a last 
stimulating chapter. There is a good bibliography of Duhem and his critics, to which should be added the 


rather important books by A. Rey and P. Humbert mentioned on pp. 15 and 8 respectively. 
P. LE CorRBEILLER 




















